import numpy as np
import matplotlib.pyplot as plt


def tabla_diferencias_divididas(x, y):
    """
    Construye la tabla de diferencias divididas de Newton.
    """
    x = np.array(x, dtype=float)
    y = np.array(y, dtype=float)

    n = len(x)
    tabla = np.zeros((n, n))
    tabla[:, 0] = y

    for j in range(1, n):
        for i in range(n - j):
            tabla[i, j] = (tabla[i + 1, j - 1] - tabla[i, j - 1]) / (x[i + j] - x[i])

    return tabla


def evaluar_newton(x, coef, valor):
    """
    Evalúa el polinomio de Newton en un punto valor.
    """
    n = len(coef)
    resultado = coef[n - 1]

    for k in range(n - 2, -1, -1):
        resultado = resultado * (valor - x[k]) + coef[k]

    return resultado


def imprimir_tabla(x, tabla):
    """
    Imprime la tabla triangular de diferencias divididas.
    """
    n = len(x)

    print("\nTabla de diferencias divididas:\n")
    print(f"{'x_i':>10}", end="")
    for j in range(n):
        print(f"{'Orden ' + str(j):>15}", end="")
    print()

    for i in range(n):
        print(f"{x[i]:>10.4f}", end="")
        for j in range(n - i):
            print(f"{tabla[i, j]:>15.6f}", end="")
        print()


def polinomio_newton_texto(x, coef):
    """
    Escribe el polinomio de Newton en forma textual.
    """
    expresion = f"{coef[0]:.6f}"

    for k in range(1, len(coef)):
        termino = f" + ({coef[k]:.6f})"
        for j in range(k):
            termino += f"(x - {x[j]:.6f})"
        expresion += termino

    return expresion


# -------------------------------
# Datos del problema
# -------------------------------

x_datos = [0, 1, 2]
y_datos = [1, 3, 7]

# -------------------------------
# Cálculo del polinomio
# -------------------------------

tabla = tabla_diferencias_divididas(x_datos, y_datos)
coeficientes = tabla[0, :]

imprimir_tabla(x_datos, tabla)

print("\nCoeficientes del polinomio de Newton:")
print(coeficientes)

print("\nPolinomio de Newton:")
print("P(x) =", polinomio_newton_texto(x_datos, coeficientes))

# -------------------------------
# Evaluación del polinomio
# -------------------------------

x_eval = 1.5
y_eval = evaluar_newton(x_datos, coeficientes, x_eval)

print(f"\nEvaluación:")
print(f"P({x_eval}) = {y_eval:.6f}")

# -------------------------------
# Gráfica
# -------------------------------

x_grafica = np.linspace(min(x_datos), max(x_datos), 200)
y_grafica = [evaluar_newton(x_datos, coeficientes, xi) for xi in x_grafica]

plt.figure(figsize=(7, 5))
plt.plot(x_grafica, y_grafica, label="Polinomio interpolador")
plt.scatter(x_datos, y_datos, label="Datos", zorder=3)
plt.xlabel("x")
plt.ylabel("P(x)")
plt.title("Interpolación polinomial de Newton")
plt.grid(True)
plt.legend()
plt.show()