# polynomial interpolation according to Neville's algorithm,
# i.e.  "Numerical Recipes" polint() on page 109 translated into Python.
# polint() takes two lists, xa, ya, size n, and input x, and pass out
# two lists of size 1 to emulate C for y and dy.  y is interpolated value
# and dy is its error. Unlike the C code, here list index starts from 0.
#
import math

def polint(xa, ya, n, x, y, dy):
    # initialize c and d to be equal to ya
 
    c = ya.copy()
    d = ya.copy()
    # then find an index ns which is closest to the value x in xa()
    ns = 0
    dif = math.fabs(x-xa[0])
    for i in range(1,n):
        dift = math.fabs(x-xa[i])
        if (dift < dif):
            ns = i
            dif = dift
    y[0] = ya[ns]
    ns -= 1

    # do double loop over column m and row i on the triangular Tableau.
    # we don't try to catch the dividing by 0 error, let Python itself do it
    for m in range(1,n):
        for i in range(0,n-m):
            ho = xa[i]-x
            hp = xa[i+m]-x
            w = c[i+1]-d[i]
            den = ho - hp
            den = w/den
            d[i] = hp*den
            c[i] = ho*den
        # end for i
        # tricking coding here in C
        if(2*(ns+1) < (n-m)):
            dy[0] = c[ns+1]
        else:
            dy[0] = d[ns]
            ns -= 1
        # end if else
        y[0] += dy[0]
    # end for m
    print(c)
    print(d)
    
# end def polint(...)

# a test run
xa = [0,1,2,4]
ya = [1,2,3,0]
x = 3.0
y = [0]
dy = [0]
n = 4
polint(xa,ya,n,x,y,dy)
print("y=",y, "dy=",dy)