Newton Forward Divided Difference Formula

The Newton polynomial can be expressed in a simplified form when $ x_0,x_1,...,x_k $ are arranged consecutively with equal spacing. Introducing the notation $h=x_{i+1} - x_i$ for each $ i = 0, 1,...,k-1$ and $ x = x_0+sh$, the difference $x-x_i$ can be written as $(s-i)h$. So the Newton polynomial becomes $$\sum_{i=0}^k = {s \choose i}i!h^i[y_0,...,y_i]. $$This is called the Newton forward divided difference formula.

Let's check this method for the next function: $$f(x*3^x)$$ with $\varepsilon = 0.0001$

import numpy as np 
import math
x = np.array([-1.0, -0.5, 0, 0.5, 1.0], float)
y = np.array([-0.3333333333333333,-0.28867513459481287, 0.0, 0.8660254037844386, 3.0], float)
eps = 0.0001
n = len(x)

for i in range(n):
    a = float(x[i]*(3**x[i]))
    print("f(x",i,") =",a)

def coef(x, y): 
    x.astype(float) 
    y.astype(float) 
    n = len(x) 
    a = [] 
    for i in range(n): 
        a.append(y[i]) 
    for j in range(1, n): 
        for i in range(n-1, j-1, -1): 
            a[i] = float(a[i]-a[i-1])
    return np.array(a)  
print(coef(x, y))
r = float(-0.25)
t = float((r - x[0])/(x[1]-x[0]))

def tsum(t, n):
    sum = float(1.0)
    for i in range(0, n, 1):
        sum = float(sum*(t-i))
    return float(sum)
    
def Eval(a, x, t): 
    x.astype(float) 
    n = len(a) 
    temp = a[0]
    count = 1
    for i in range(1, n, 1):
        temp += (tsum(t, i-1)*a[i])/math.factorial(i)
        count = count + 1
    print("Count of parts in the interpolation formula:",count)
    return temp
result = Eval(coef(x,y),x,t)

print("Value of functions f(x) in",r,"=",result)
f(x 0 ) = -0.3333333333333333
f(x 1 ) = -0.28867513459481287
f(x 2 ) = 0.0
f(x 3 ) = 0.8660254037844386
f(x 4 ) = 3.0
[-0.33333333  0.0446582   0.24401694  0.33333333  0.35726559]
Count of parts in the interpolation formula: 5
Value of functions f(x) in -0.25 = -0.06957804087824197