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

Question

我将在此处和下面的上下文中进行简要说明。我在 Python 和 Fortran 中评估双变量多项式(2 个变量的多项式)并得到不同的结果。我的测试用例 - 4.23e-3 - 的相对误差足够大，不会因为精度差异而明显。以下代码片段使用相当原始的类型和相同的算法来尝试使事物尽可能具有可比性。关于差异的任何线索？我尝试过改变精度(在 Fortran 中使用 selected_real_kind 和在 Python 中使用 numpy.float128)，Fortran 的编译(特别是优化级别)和算法(霍纳的方法，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

---

语言:

! 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 实现，主要是作为练习，但也为了增加一些灵 active 。我正在将我的结果与 Fortran 进行基准测试，并且已经能够在大约 1e-10 内获得几乎所有内容，但是这个超出了我的掌握范围。其他函数也复杂得多，使得简单多项式的分歧令人费解。

特定的系数和测试变量来自这个库。实际的多项式实际上在 (x,y) 中有 (7,6) 次，所以这里没有包括更多的系数。算法直接取自Fortran，所以如果有错，我应该联系原开发者。一般函数也可以计算导数，这就是为什么这个实现可能不是最优的部分原因——我知道我仍然应该只写一个 Horner 的方法版本，但这并没有改变差异。我只在计算大 y 值的导数时才注意到这些错误，但错误确实持续存在于这个更简单的设置中。

Answer 1

应更正 Fortran 代码中的两件事，以使结果在 Python 和 Fortran 版本之间匹配。

1.正如您所做的那样，将特定的 double 类型声明为:

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

然后您应该通过附加类型指示符来定义一个变量:

real(dbl) :: z
z = 1.0_dbl

例如，在 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