Accuracy and polynomial evaluation in Fortran and Python

Question

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.

Answer 1

Two things in the Fortran code should be corrected to get the results to match between the Python and Fortran versions.

1. As you have done, declare a specific double-precision kind as:

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

You should then define a variable by appending the kind designator as:

real(dbl) :: z
z = 1.0_dbl

This is discussed, for example, at fortran90.org gotchas. The syntax may be inconvenient, but hey, I don't make the rules.

2. The Fortran do-loop iteration is controlled by nx and ny. You intend to access each element of coeffs, but your indexing is cutting the iteration short. Change nx-1 and ny-1 to nx and ny, respectively. Even better, use the Fortran intrinsic ubound to determine the extent along the desired dimension, such as:

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

The updated code shown below corrects these issues and prints the result, which matches that produced by your python code.

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