#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jun  6 14:10:02 2023

@author: Christoph Matern
"""

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import barycentric_interpolate
from copy import copy

# %%


def _newton_divdiff(x, y):
    """
    siehe  https://stackoverflow.com/questions/14823891/newton-s-interpolating-polynomial-python
    x: list or np array contanining x data points
    y: list or np array contanining y data points
    """
    m = len(x)
    dd = copy(y)

    for k in range(1, m):
        dd[k:m] = (dd[k:m] - dd[k - 1])/(x[k:m] - x[k - 1])

    return dd

# %%


def newton_interpolate(x, y, X):
    """
    x: data points at x
    y: data points at y
    X: evaluation point(s)
    """
    d = _newton_divdiff(x, y)
    n = len(x) - 1  # Grad des Polynoms
    p = d[n]

    for k in range(1, n + 1):
        p = d[n - k] + (X - x[n - k])*p

    return p

# %%


def polynome(m, L):

    # Stützstellen als rationale Zahlen
    x_s = np.array([sp.S(L)/m*j*2-L for j in range(m+1)])
    yrand = np.random.randint(-123, 123, len(x_s))
    y_s = np.array(sp.S(yrand)/123)

    # Stützstellen als Gleitkommazahl mit doppelter Genauigkeit
    x_n = np.float64(x_s)
    y_n = np.float64(y_s)

    def newton_polynom_n(X):
        return newton_interpolate(x_n, y_n, X)

    x = sp.symbols('x')
    # Die Hornerform ist etwas übersichtlicher
    p = sp.horner(sp.interpolating_poly(len(x_s), x, X=x_s, Y=y_s))
    if m < 13:
        print(p)

    def p_exakt(X):
        return np.float64(sp.N(p.subs(x, X), 24))

    def bary_polynom_n(X):
        pass

    def vandermonde_polynom_n(X):
        pass

    return newton_polynom_n, bary_polynom_n, vandermonde_polynom_n, p_exakt

# %%


m = 12
L = 2
p_newton, p_bary, p_vandermonde, p_exakt = polynome(m, L)

err_n = []
nref = 123
xs = [sp.Rational(L, nref)*j for j in range(-nref, nref+1)]
xn = np.float64(xs)

for xn_, xs_ in zip(xn, xs):
    err_n.append(p_newton(xn_) - p_exakt(xs_))

plt.plot(xn, err_n, label='Newton')
plt.xlabel("x")
plt.ylabel('Fehler')
plt.legend()