Экспонирование в Фортран / Гфортран с высокой точностью - PullRequest
Новая
       52

Экспонирование в Фортран / Гфортран с высокой точностью

2 голосов
/ 10 июля 2019

Как gfortran обрабатывает возведение в степень с целым числом против действительного?Я всегда предполагал, что это было то же самое, но рассмотрим пример: 

program main

implicit none

integer, parameter :: dp = selected_real_kind(33,4931)

real(kind=dp) :: x = 82.4754500815524510_dp

print *, x
print *, x**4
print *, x**4.0_dp

end program main

Компиляция с gfortran дает 

82.4754500815524510000000000000000003      
46269923.0191143410452125643548442147      
46269923.0191143410452125643548442211

Теперь ясно, что эти числа почти совпадают - но если gfortran обрабатывает целые числа и выдает реал длявозведение в степень так же, как я ожидал бы, что они будут идентичны.Что дает?

Ответы [ 2 ]

5 голосов
/ 10 июля 2019

Расширение вашей программы немного показывает, что происходит: 

ijb@ianbushdesktop ~/work/stack $ cat exp.f90
program main

implicit none

integer, parameter :: dp = selected_real_kind(33,4931)

real(kind=dp) :: x = 82.4754500815524510_dp

print *, x
print *, x**4
print *,(x*x)*(x*x)
print *,Nearest((x*x)*(x*x),+1.0_dp)
print *, x**4.0_dp

end program main

Компиляция и запуск дает: 

ijb@ianbushdesktop ~/work/stack $ gfortran --version
GNU Fortran (GCC) 7.4.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ijb@ianbushdesktop ~/work/stack $ gfortran exp.f90 
ijb@ianbushdesktop ~/work/stack $ ./a.out
   82.4754500815524510000000000000000003      
   46269923.0191143410452125643548442147      
   46269923.0191143410452125643548442147      
   46269923.0191143410452125643548442211      
   46269923.0191143410452125643548442211      
ijb@ianbushdesktop ~/work/stack $

Таким образом

  1. Этопохоже, что компилятор достаточно умен, чтобы понять, что целочисленные полномочия могут быть получены умножением, что намного быстрее, чем общая функция возведения в степень

  2. Использование общей функции возведения в степень отличается только на 1 битиз умножения ответа.Поскольку мы не можем сказать per se , что является более точным, мы должны принять оба одинаково точных

Таким образом, в заключение, компилятор использует простые умножения, где это возможно, вместополномасштабная процедура возведения в степень, но даже если она использует более дорогой маршрут, она дает тот же ответ, для тщательно продуманного значения слова «то же самое»

1 голос
/ 10 июля 2019

Для полноты, возведение в степень (точной) дроби приводит к 

46269923.019114341045212564354844220930226304938209797955723262974801
46269923.0191143410452125643548442211 is the nearest
46269923.0191143410452125643548442147

Возведение в степень ближайшей четвертой точности с плавающей точкой (https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format) приводит к 

46269923.01911434104521256435484422112355946434320203876355837024902939447201788425445556640625
46269923.0191143410452125643548442211 does fit

Похоже, функция возведения в степень округлилась до ближайшего числа с плавающей точкой. Я не знаю, является ли это удачей или если соответствующая математическая библиотека гарантировала правильное округление.

Целочисленное экспонента: 

x2 = x*x
x4 = x2*x2

накапливает ошибку двух округлений, так что это может быть 1 отключение, это неудивительно.

...