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