floating-point - 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. 看起来编译器足够聪明，可以计算出整数幂可以通过乘法来完成，这比一般的求幂函数快得多

2. 使用一般求幂函数与乘法答案只有 1 位不同。鉴于我们不能说本身哪个更准确，我们必须接受两者同样准确

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