Fortran和Python中的精度和多项式评估

Question

我将在这里和下面的上下文中给出一个快速说明。我正在Python和Fortran中评估双变量多项式（2个变量的多项式）并得到不同的结果。我的测试用例 - 4.23e-3的相对误差足够大，不会因为精度差异而显着。下面的代码片段使用相当原始的类型和相同的算法来尝试尽可能地比较。任何关于差异的线索？我尝试了各种精度（包括Fortran中的selected_real_kind和Python中的numpy.float128），Fortran编译（特别是优化级别）和算法（Horner的方法，numpy评估）。任何关于差异的线索？代码版本中的错误？我见过Precision discrepancy between Fortran and Python (sin function)，但没有机会完全使用不同的编译器进行测试。

---

的Python：

#!/usr/bin/env python 
""" polytest.py 
Test calculation of a bivariate polynomial. 
""" 

# Define polynomial coefficients 
coeffs = (
    (101.34274313967416, 100015.695367145, -2544.5765420363), 
    (5.9057834791235253,-270.983805184062,1455.0364540468), 
    (-12357.785933039,1455.0364540468,-756.558385769359), 
    (736.741204151612,-672.50778314507,499.360390819152) 
) 
nx = len(coeffs) 
ny = len(coeffs[0]) 

# Values of variables 
x0 = 0.0002500000000011937 
y0 = -0.0010071334522899211 

# Calculate polynomial by looping over powers of x, y 
z = 0. 
xj = 1. 
for j in range(nx): 
    yk = 1. 
    for k in range(ny): 
     curr = coeffs[j][k] * xj * yk 
     z += curr 
     yk *= y0 
    xj *= x0 

print(z) # 0.611782174444 

---

的Fortran：

! polytest.F90 
! Test calculation of a bivariate polynomial. 

program main 

    implicit none 
    integer, parameter :: dbl = kind(1.d0) 
    integer, parameter :: nx = 3, ny = 2 
    real(dbl), parameter :: x0 = 0.0002500000000011937, & 
          y0 = -0.0010071334522899211 
    real(dbl), dimension(0:nx,0:ny) :: coeffs 
    real(dbl) :: z, xj, yk, curr 
    integer :: j, k 

    ! Define polynomial coefficients 
    coeffs(0,0) = 101.34274313967416d0 
    coeffs(0,1) = 100015.695367145d0 
    coeffs(0,2) = -2544.5765420363d0 
    coeffs(1,0) = 5.9057834791235253d0 
    coeffs(1,1) = -270.983805184062d0 
    coeffs(1,2) = 1455.0364540468d0 
    coeffs(2,0) = -12357.785933039d0 
    coeffs(2,1) = 1455.0364540468d0 
    coeffs(2,2) = -756.558385769359d0 
    coeffs(3,0) = 736.741204151612d0 
    coeffs(3,1) = -672.50778314507d0 
    coeffs(3,2) = 499.360390819152d0 

    ! Calculate polynomial by looping over powers of x, y 
    z = 0d0 
    xj = 1d0 
    do j = 0, nx-1 
    yk = 1d0 
    do k = 0, ny-1 
     curr = coeffs(j,k) * xj * yk 
     z = z + curr 
     yk = yk * y0 
    enddo 
    xj = xj * x0 
    enddo 

    ! Print result 
    WRITE(*,*) z ! 0.61436839888538231 

end program 

编译时：gfortran -O0 -o polytest.o polytest.F90

---

上下文：我正在编写一个现有Fortran库的纯Python实现，主要是作为练习，但也增加了一些灵活性。我将我的结果与Fortran进行了比较，并且能够在大约1e-10的范围内获得几乎所有的结果，但是这一点在我的掌握之外。其他函数也更加复杂，使得简单多项式的分歧令人困惑。

特定的系数和测试变量来自这个库。实际的多项式实际上在（x，y）中有度（7,6），因此这里没有包含更多的系数。该算法直接来自Fortran，所以如果它是错误的，我应该联系原始开发人员。一般函数也可以计算衍生物，这是这种实现可能不是最优的原因之一 - 我知道我应该只写Horner的方法版本，但这并没有改变这种差异。在计算y的大数值时，我只注意到这些错误，但错误仍然存​​在于这个更简单的设置中。

Answer 1

Fortran代码中的两件事情应该得到纠正以获得匹配Python和Fortran版本的结果。

正如你所做的那样，宣布具体的双精度种为：

real(dbl) :: z 
z = 1.0_dbl 

此：

integer, parameter :: dbl = kind(0.d0) 

然后，您应该通过附加一种代号为定义一个变量例如在fortran90.org gotchas中讨论。语法可能不方便，但是，嘿，我没有制定规则。

2. Fortran do-loop迭代由nx和ny控制。您打算访问coeffs的每个元素，但是您的索引将缩短迭代次数。分别将nx-1和ny-1更改为nx和ny。更好的是，使用Fortran固有ubound程度来决定沿着期望的尺寸，例如：

do j = 0, ubound(coeffs, dim=1) 

下面显示的更新的代码纠正这些问题，并输出结果，它匹配，通过Python代码产生的。

program main 
    implicit none 
    integer, parameter :: dbl = kind(1.d0) 
    integer, parameter :: nx = 3, ny = 2 
    real(dbl), parameter :: x0 = 0.0002500000000011937_dbl, & 
          y0 = -0.0010071334522899211_dbl 
    real(dbl), dimension(0:nx,0:ny) :: coeffs 
    real(dbl) :: z, xj, yk, curr 
    integer :: j, k 

    ! Define polynomial coefficients 
    coeffs(0,0) = 101.34274313967416_dbl 
    coeffs(0,1) = 100015.695367145_dbl 
    coeffs(0,2) = -2544.5765420363_dbl 
    coeffs(1,0) = 5.9057834791235253_dbl 
    coeffs(1,1) = -270.983805184062_dbl 
    coeffs(1,2) = 1455.0364540468_dbl 
    coeffs(2,0) = -12357.785933039_dbl 
    coeffs(2,1) = 1455.0364540468_dbl 
    coeffs(2,2) = -756.558385769359_dbl 
    coeffs(3,0) = 736.741204151612_dbl 
    coeffs(3,1) = -672.50778314507_dbl 
    coeffs(3,2) = 499.360390819152_dbl 

    ! Calculate polynomial by looping over powers of x, y 
    z = 0.0_dbl 
    xj = 1.0_dbl 
    do j = 0, ubound(coeffs, dim=1) 
     yk = 1.0_dbl 
     do k = 0, ubound(coeffs, dim=2) 
      print "(a,i0,a,i0,a)", "COEFF(",j,",",k,")=" 
      print *, coeffs(j,k) 
      curr = coeffs(j,k) * xj * yk 
      z = z + curr 
      yk = yk * y0 
     enddo 
     xj = xj * x0 
    enddo 

    ! Print result 
    WRITE(*,*) z ! Result: 0.611782174443735 
end program