# newtonInterp.py
# Andrew Davison, Sept. 2025, ad@coe.psu.ac.th
'''
Polynomial Interpolation:

Find a polynomial P_n(x) of degree n that passes exactly 
through every data point (x_i, y_i).
Condition: number of points = degree + 1 for a unique polynomial.

Using Newton's Divided Difference Method 

Build the polynomial incrementally in Newton's form: 

P(x) = a0 + a1(x - x0) + a2(x - x0)(x - x1) + ...
       + an(x - x0)(x - x1)...(x - xn_1) 

where the ak are computed using divided differences. 

Compute coefficients with a triangular divided-difference table. 
Evaluate using nested multiplication.
'''


import math
import matplotlib.pyplot as plt


def evalPoly(a, xData, x):
  # Evaluate Newton's polynomial using nested multiplication
  n = len(xData) - 1
  poly = a[n]
  for k in range(1, n + 1):
    poly = a[n - k] + (x - xData[n - k]) * poly
  return poly


def coefs(xData, yData):
  # Compute divided difference coefficients
  m = len(xData)
  a = yData[:]  # make a copy
  for k in range(1, m):
    for j in range(m-1, k-1, -1):
      a[j] = (a[j] - a[j-1]) / (xData[j] - xData[j-k])
  return a


def polyEquation(a, xData):
  # Build human-readable polynomial equation string
  terms = [f"{a[0]:.5f}"]
  for k in range(1, len(a)):
    factors = "".join([f"(x-{xData[j]:.2f})" for j in range(k)])
    terms.append(f"{a[k]:+0.3f}{factors}")
  return " ".join(terms)


# Example data
xData = [0.15, 2.3, 3.15, 4.85, 6.25, 7.95]
yData = [4.79867, 4.49013, 4.2243, 3.47313, 2.66674, 1.51909]

a = coefs(xData, yData)

# Print polynomial equation to stdout
print("Newton Interpolation Polynomial:")
print("P(x) =", polyEquation(a, xData))
print()

# Build points for plotting
xs = []
ysInterp = []
ysExact = []
x = 0.0
while x <= 8.0:
  xs.append(x)
  ysInterp.append(evalPoly(a, xData, x))
  ysExact.append(4.8 * math.cos(math.pi * x / 20.0))
  x += 0.1  # smoother curve

# Plot
plt.figure(figsize=(8, 5))
plt.plot(xs, ysInterp, "r-", label="Newton Polynomial Interpolation")
plt.plot(xs, ysExact, "b--", label="Exact Function $4.8\\cos(\\pi x/20)$")
plt.scatter(xData, yData, c="k", marker="o", label="Data Points")

plt.xlabel("x")
plt.ylabel("y")
plt.title("Newton Interpolation Polynomial vs Exact Function")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()