import numpy as np
import matplotlib.pyplot as plt

px = py = []
with open("data8//data15.txt") as file:
    for line in file.readlines():
        lines = line.split()
        px = np.append(px, float(lines[0].strip()))
        py = np.append(py, float(lines[1].strip()))

n = len(px)
f = {}
for i in range(0, n):
    f[px[i]] = py[i]


def fun1(x2, x1):
    return (f.get(x2) - f.get(x1)) / (x2 - x1)


def fun2(x3, x2, x1):
    return (fun1(x3, x2) - fun1(x2, x1)) / (x3 - x1)


def fun3(x4, x3, x2, x1):
    return (fun2(x4, x3, x2) - fun2(x3, x2, x1)) / (x4 - x1)


def fun4(x5, x4, x3, x2, x1):
    return (fun3(x5, x4, x3, x2) - fun3(x4, x3, x2, x1)) / (x5 - x1)


def fun5(x6, x5, x4, x3, x2, x1):
    return (fun4(x6, x5, x4, x3, x2) - fun4(x5, x4, x3, x2, x1)) / (x6 - x1)


def fun6(x7, x6, x5, x4, x3, x2, x1):
    return (fun5(x7, x6, x5, x4, x3, x2) - fun5(x6, x5, x4, x3, x2, x1)) / (x7 - x1)


def fun7(x8, x7, x6, x5, x4, x3, x2, x1):
    return (fun6(x8, x7, x6, x5, x4, x3, x2) - fun6(x7, x6, x5, x4, x3, x2, x1)) / (x8 - x1)


def fun8(x9, x8, x7, x6, x5, x4, x3, x2, x1):
    return (fun7(x9, x8, x7, x6, x5, x4, x3, x2) - fun7(x8, x7, x6, x5, x4, x3, x2, x1)) / (x9 - x1)


def fun9(x10, x9, x8, x7, x6, x5, x4, x3, x2, x1):
    return (fun8(x10, x9, x8, x7, x6, x5, x4, x3, x2) - fun8(x9, x8, x7, x6, x5, x4, x3, x2, x1)) / (x10 - x1)


def fun10(x11, x10, x9, x8, x7, x6, x5, x4, x3, x2, x1):
    return (fun9(x11, x10, x9, x8, x7, x6, x5, x4, x3, x2) - fun9(x10, x9, x8, x7, x6, x5, x4, x3, x2, x1)) / (x11 - x1)

b = [1]
b.append(f.get(px[0]))
b.append(fun1(px[1], px[0]))
b.append(fun2(px[2], px[1], px[0]))
b.append(fun3(px[3], px[2], px[1], px[0]))
b.append(fun4(px[4], px[3], px[2], px[1], px[0]))
b.append(fun5(px[5], px[4], px[3], px[2], px[1], px[0]))
b.append(fun6(px[6], px[5], px[4], px[3], px[2], px[1], px[0]))
b.append(fun7(px[7], px[6], px[5], px[4], px[3], px[2], px[1], px[0]))
b.append(fun8(px[8], px[7], px[6], px[5], px[4], px[3], px[2], px[1], px[0]))
b.append(fun9(px[9], px[8], px[7], px[6], px[5], px[4], px[3], px[2], px[1], px[0]))
b.append(fun10(px[10], px[9], px[8], px[7], px[6], px[5], px[4], px[3], px[2], px[1], px[0]))

def newton(x):
    factor = b[1]
    out = factor
    for i in range(2, n+1):
        factor = factor*(b[i]/b[i-1])*(x-px[i-2])
        out += factor
    return out

i = np.linspace(min(px), max(px), 200)
plt.plot(px, py, 'r.', i, newton(i))
plt.show()