#!/usr/bin/env python3
"""Compare degree-n Runge interpolation on equispaced and Lobatto nodes.

Run with Python 3 and NumPy. The reported error is a maximum on a fixed
20,001-point grid, not a certified continuous supremum. Floating-point
summation and the platform math library can change the final digits.
"""
import json
import math
import numpy as np


def barycentric(nodes, values, weights, evaluation_points):
    """Second barycentric formula, with exact-node values restored."""
    differences = evaluation_points[:, None] - nodes[None, :]
    exact = differences == 0
    safe = np.where(exact, 1.0, differences)
    terms = weights[None, :] / safe
    result = (terms @ values) / terms.sum(axis=1)
    rows, columns = np.where(exact)
    result[rows] = values[columns]
    return result


def main():
    grid = np.linspace(-1.0, 1.0, 20001)
    truth = 1.0 / (1.0 + 25.0 * grid**2)
    rows = []
    for degree in (10, 20, 30):
        equispaced = np.linspace(-1.0, 1.0, degree + 1)
        equal_weights = np.array([
            (-1.0)**j * math.comb(degree, j)
            for j in range(degree + 1)
        ])
        lobatto = np.cos(np.arange(degree + 1) * np.pi / degree)
        lobatto_weights = (-1.0)**np.arange(degree + 1)
        lobatto_weights[[0, -1]] *= 0.5
        errors = []
        for nodes, weights in ((equispaced, equal_weights),
                               (lobatto, lobatto_weights)):
            values = 1.0 / (1.0 + 25.0 * nodes * nodes)
            approximation = barycentric(nodes, values, weights, grid)
            errors.append(float(np.max(np.abs(approximation - truth))))
        rows.append({"degree": degree, "equispaced_error": errors[0],
                     "lobatto_error": errors[1], "ratio": errors[0] / errors[1]})
    print(json.dumps({"numpy_version": np.__version__, "grid_points": len(grid),
                      "rows": rows}, indent=2))


if __name__ == "__main__":
    main()