from math import *
import numpy as np
import matplotlib.pyplot as plt


def f(x):
    y = 1 / (1 + 25 * x**2)
    return y


def diffdiv(x, y):
    n = len(x) - 1
    D = np.zeros((n+1, n+1))
    D[:,0] = y
    for j in range (1, n+1):
        for i in range (n-j+1):
            D[i,j] = (D[i+1,j-1] - D[i,j-1]) / (x[i+j] - x[i])
    return D

def fastexp(x, D, t):
    n = len(x) - 1
    y = D[0, n] * np.ones(len(t))
    for k in range (n-1, -1, -1):
        y = (t-x[k]) * y + D[0, k]
    return y


x = np.linspace(-1, 1, 400)
y = f(x)
plt.plot(x, y, label = "f")

for n in range (1, 10):
    X = np.linspace(-1, 1, n+1)
    Y = f(X)
    D = diffdiv(X, Y)
    y = fastexp(X, D, x)
    plt.plot(x, y, '--', label = "n = " + str(n))
    plt.plot(X, Y, 'Xr')
plt.legend(loc = "best")
plt.title("Polynômes d'interpolation de degré n")
plt.show()