I'll give a quick description here and the context below. I'm evaluating bivariate polynomials (polynomials of 2 variables) in both Python and Fortran and getting different results. The relative error for my test case - 4.23e-3 - is large enough to not be obviously due to precision differences. The following code snippets use fairly primitive types and the same algorithm to try to make things as comparable as possible. Any clues as to the discrepancy? I have tried varying the precision (both in Fortran with selected_real_kind
and in Python with numpy.float128
), compilation of the Fortran (specifically, the optimization level), and the algorithm (Horner's method, numpy
evaluation). Any clues as to the discrepancy? Errors in either version of the code? I've seen Precision discrepancy between Fortran and Python (sin function) but haven't had a chance to test it with different compilers entirely.
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
Compiled with: gfortran -O0 -o polytest.o polytest.F90
Context: I'm writing a pure-Python implementation of an existing Fortran library, primarily as an exercise but also to add some flexibility. I'm benchmarking my results against the Fortran and have been able to get almost everything within about 1e-10, but this one is outside of my grasp. The other functions are also much more complex, making the disagreement for simple polynomials baffling.
The particular coefficients and test variables come from this library. The actual polynomial actually has degree (7,6) in (x,y), so there are many more coefficients not included here. The algorithm is taken directly from the Fortran, so if it's wrong, I should contact the original developers. The general functions can also calculate the derivatives, which is part of why this implementation is maybe sub-optimal -- I'm aware that I should still just write a Horner's method version, but that didn't change the discrepancy. I only noticed these errors when calculating the derivatives at large values of y, but the error does persist into this simpler setting.
j
is going to go from0 -> 2 (nx - 1)
. In the python code,j
is going from0 -> 3 (len(coeffs) - 1)
– mgilson