Newton Polynomial in Python

I wrote code for this in Delphi. This time I want to rewrite it in Python based on this wiki.

I use this set of data point
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)

and I use xc=3 for the test data.

It's obvious that these sets of data points have quadratic form and f(xc) must have value of 9.

The heart of code lay on this

def n(j,xc,x):
    n = 1
    for i in arange(j):
        n *= (xc-x[i])

    return n
def a(j,l,x,y):
    if j==0:
        return y[0]
    elif j-l==1 :
        return (y[j]-y[l])/(x[j]-x[l])
    else:
        return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])
        

def N(xc,x,y):
    N = 0
    for j in range(len(x)):
        N += a(j,0,x,y)*n(j,xc,x)
    
    return N

Look at the function a(j,l,x,y), that's recursive function to obtain Newton divided difference value.

the whole code is this

from pylab import *

def n(j,xc,x):
    n = 1
    for i in arange(j):
        n *= (xc-x[i])

    return n
def a(j,l,x,y):
    if j==0:
        return y[0]
    elif j-l==1 :
        return (y[j]-y[l])/(x[j]-x[l])
    else:
        return (a(j,l+1,x,y)-a(j-1,l,x,y))/(x[j]-x[l])
        

def N(xc,x,y):
    N = 0
    for j in range(len(x)):
        N += a(j,0,x,y)*n(j,xc,x)
    
    return N

x = []
y = []
#initial value
x.append(0)
x.append(1)
x.append(2)
x.append(4)
x.append(5)

y.append(0)
y.append(1)
y.append(4)
y.append(16)
y.append(25)

#for testing
xc = 3

yc = N(xc,x,y)

print ''
print xc, yc
#plot
t = linspace(-7,7,100)
u = N(t,x,y)
plot(t,u)
grid(True)
show()

.


here's the graphics

with another sets of data points, I have another result
(0,0)
(1,1)
(2,4)
(4,16)
(5,25)