Fortran/gfortran 中的求幂到高精度

Question

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 以相同的方式处理整数和实数以求幂，我希望它们是相同的。给出了什么？

Answer 1

稍微扩展您的程序可以看出发生了什么：

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位不同。鉴于我们不能说本身哪个更准确，我们必须接受两者同样准确

因此，编译器在可能的情况下使用简单的乘法而不是完整的求幂例程，但即使它必须使用更昂贵的路径，它也会给出相同的答案，因为仔细考虑了这个词的含义"same"

Answer 2

为了完整起见，对（精确）分数取幂得到

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 ulp，这不足为奇。

Fortran/gfortran 中的求幂到高精度

Exponentiation in Fortran/gfortran to high precision

floating-point

fortran

gfortran

exponentiation