Precision and polynomial estimation in Fortran and Python

I'll provide a short description here and context below. I evaluate 2D polynomials (polynomials in 2 variables) in both Python and Fortran and getting different results. The relative error for my test - 4.23e-3 - is large enough not to be obvious due to differences in accuracy. The following code snippets use pretty primitive types and the same algorithm to try and make things as comparable as possible. Any clues regarding the discrepancy? I have tried changing the precision (in both Fortran c selected_real_kind

and Python c numpy.float128

), Fortran compilation (specifically the optimization level), and the algorithm (Horner's method numpy

). Any clues regarding the discrepancy? Errors in any version of the code? I saw the exact mismatch between Fortran and Python (sin function)but he hasn't had the opportunity to fully test it with different compilers.

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 to add some flexibility. I compare my results with Fortran and get pretty much everything in the 1-10 range, but that's not in my hands. Other functions are also much more complex, which leads to inconsistency of simple polynomials.

Specific coefficients and test variables come from this library. The actual polynomial actually has degree (7.6) in (x, y), so there aren't many coefficients here. The algorithm is taken directly from Fortran, so if this is wrong I should contact the original developers. Generic functions can also compute derivatives, which is part of why this implementation might not be optimal - I know I still need to write a version of Horner's method, but that hasn't changed the discrepancy. I only noticed these errors when calculating derivatives at large y values, but the error persists in this simpler setting.