【发布时间】：2017-08-03 12:25:21
【问题描述】：
我将在此处和下面的上下文中进行简要说明。我正在评估 Python 和 Fortran 中的二元多项式（2 个变量的多项式）并得到不同的结果。我的测试用例的相对误差 - 4.23e-3 - 足够大，不会因为精度差异而明显。以下代码 sn-ps 使用相当原始的类型和相同的算法来尝试使事物尽可能具有可比性。关于差异的任何线索？我尝试改变精度（在 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，如有不对请联系原开发者。通用函数还可以计算导数，这就是为什么这个实现可能不是最佳的部分原因——我知道我仍然应该只写一个霍纳的方法版本，但这并没有改变差异。我只在计算大 y 值的导数时才注意到这些错误，但在这个更简单的设置中，错误确实存在。
【问题讨论】：
-
如果它们都是同一个算法，你应该能够比较两个实现之间的中间值，看看差异是如何发展的，如果有“罪魁祸首”，找出它。跨度>
-
我已经很久没有看过 Fortran 代码了，但看起来你的循环边界不同。例如，在 fortran 代码中，
j将来自0 -> 2 (nx - 1)。在 python 代码中，j来自0 -> 3 (len(coeffs) - 1) -
同样在 Fortran 中，您只需将参数 x0 和 y0 设置为单精度。此外，常量的规范是非常老式（f77）风格和更现代的使用类型值的奇怪组合。
-
另外，我认为@mgilson 发现了错误。
标签： python fortran precision floating-accuracy polynomial-math