Interpolation

Why Interpolation Matters

A calibration table records a sensor's output at several known inputs. A new input falls between two recorded positions. Interpolation constructs a function that passes through the recorded pairs so that the missing value can be estimated consistently. Animation curves, tabulated physical models, and resampling all use this pattern.

This lesson covers Lagrange interpolation, Newton divided differences, and cubic splines. Lagrange and Newton are two representations of the same global polynomial for a given set of distinct nodes. A cubic spline chooses a piecewise polynomial under additional smoothness and boundary conditions, so it can give different values between the same knots. The implementations use a natural cubic spline, whose endpoint second derivatives are zero.

You will construct each representation, trace its coefficients and evaluations, implement all three in Python, C++, and Java, and check knot residuals, polynomial agreement, and the spline equations. Interpolation fits the supplied values exactly in the mathematical model; it does not establish that those values are noise-free or that the curve predicts an unknown function accurately between them.

Core Concepts and Terminology

Let there be $n\ge2$ data pairs $(x_i,y_i)$ with distinct real nodes. A unique polynomial $P$ of degree at most $n-1$ satisfies all $n$ interpolation conditions. Lagrange's basis expresses those conditions directly, while Newton's basis organizes the same polynomial by divided differences.

A cubic spline instead uses a polynomial of degree at most three on each interval $[x_i,x_{i+1}]$. It interpolates both endpoints and has continuous first and second derivatives across interior knots. Natural endpoint conditions are $S''(x_0)=S''(x_{n-1})=0$. Other boundary conditions generally produce another spline.

Use the small data set $(0,0),(1,1),(2,0),(3,1)$. Its global polynomial and natural spline both pass through all four knots, but at $q=1/2$ their values are $P(q)=1$ and $S(q)=3/4$. Neither value follows from the knots alone without choosing the interpolation model.

All three supplied implementations require the same data format: at least two finite $[x,y]$ rows with strictly increasing nodes and a finite representable total span. This shared ordering is useful for interval lookup in the spline. Polynomial formulas themselves do not require increasing order mathematically; changing the order changes Newton's coefficients but not the exact interpolating polynomial.

How Interpolation Works, Step by Step

1. Validate and copy the data. Reject mismatched row shapes, nonfinite values, duplicate nodes, and unsorted nodes. Do not merge close nodes using an arbitrary tolerance: distinct coordinates may represent distinct measurements, and near-coincident nodes need a conditioning decision rather than silent replacement. The code also rejects unsupported arithmetic overflow or weight underflow.

2. Construct the Lagrange polynomial. Define

$$ L_i(q)=\prod_{j\ne i}\frac{q-x_j}{x_i-x_j},\qquad P(q)=\sum_{i=0}^{n-1}y_iL_i(q). $$

At a knot, all but one cardinal basis value vanish. Directly evaluating all products costs $O(n^2)$ per query. The programs instead precompute barycentric weights and use the equivalent second barycentric form

$$ P(q)=\frac{\sum_i w_i y_i/(q-x_i)}{\sum_i w_i/(q-x_i)},\qquad q\ne x_i. $$

If the query exactly equals a knot, return its recorded value before forming a denominator. A merely nearby query still uses the formula. A common nonzero scaling of all weights cancels from the ratio; the code normalizes them by their largest magnitude. During evaluation it also multiplies every term by $\min_i\lvert q-x_i\rvert$, avoiding a large reciprocal when the query is close to a knot. These changes improve arithmetic behavior without changing the exact polynomial.

3. Construct Newton divided differences. Start with $d_i^{(0)}=y_i$ and define

$$ d_i^{(k)}=\frac{d_i^{(k-1)}-d_{i-1}^{(k-1)}}{x_i-x_{i-k}},\qquad i\ge k. $$

The Newton coefficients are $c_k=d_k^{(k)}=f[x_0,\ldots,x_k]$, giving

$$ P(q)=c_0+c_1(q-x_0)+c_2(q-x_0)(q-x_1)+\cdots +c_{n-1}\prod_{j=0}^{n-2}(q-x_j). $$

The code stores only one coefficient array. At order $k$, update indices from $n-1$ down to $k$ so both numerator entries still belong to the previous order. Updating left to right would overwrite a value before it is used.

Evaluate Newton's form with a nested recurrence: start $v=c_{n-1}$, then for $i=n-2,\ldots,0$ replace $v$ by $c_i+(q-x_i)v$. This uses $O(n)$ arithmetic and avoids separately expanding monomials or multiplying every basis product.

4. Construct a natural cubic spline through its second derivatives. Let $h_i=x_{i+1}-x_i\gt0$, $s_i=(y_{i+1}-y_i)/h_i$, and $M_i=S''(x_i)$. Set $M_0=M_{n-1}=0$. Continuity of the first derivative at each interior knot yields

$$ h_{i-1}M_{i-1}+2(h_{i-1}+h_i)M_i+h_iM_{i+1} =6(s_i-s_{i-1}),\qquad 1\le i\le n-2. $$

This is a tridiagonal system for the interior second derivatives. The teaching code explicitly builds it in nested arrays and solves it using dense Gaussian elimination with partial row pivoting. It scales the matrix and right-hand side by the matrix's maximum entry magnitude, selects the largest available pivot magnitude in each column, swaps the entire matrix row and matching right-hand-side value, eliminates below, and back-substitutes.

The dimensionless pivot tolerance, default $10^{-12}$, rejects a normalized pivot at or below the threshold. It is a numerical acceptance policy, not proof that valid knots are mathematically singular. The natural-spline interior matrix is symmetric positive definite and strictly diagonally dominant for positive interval widths; a specialized tridiagonal solver can exploit that structure in linear time. The dense version is included to make the solve and pivoting explicit, with its higher cost stated below.

For each interval, convert the second derivatives to local cubic coefficients:

$$ \begin{aligned} a_i&=y_i,\\ b_i&=s_i-\frac{h_i(2M_i+M_{i+1})}{6},\\ c_i&=M_i/2,\\ d_i&=\frac{M_{i+1}-M_i}{6h_i}. \end{aligned} $$

The spline's $c_i$ here denotes a local quadratic coefficient, while Newton's $c_i$ above denotes a divided difference; they are stored separately in the code. With just two knots, there are no interior unknowns and the natural spline is the straight line between them.

5. Evaluate with an explicit domain policy. Lagrange and Newton accept finite queries inside or outside the knot range, subject to arithmetic checks. The natural-spline method rejects queries outside $[x_0,x_{n-1}]$. For an accepted spline query, binary search chooses its interval, then evaluate $((d_it+c_i)t+b_i)t+a_i$. At the final knot the search selects the final interval.

flowchart TD
    A["Validate and copy ordered knot rows"] --> B["Build normalized barycentric weights"]
    B --> C["Build Newton coefficients with descending updates"]
    C --> D{"Interior spline unknowns?"}
    D -->|Yes| E["Assemble system; pivot, eliminate, back-substitute"]
    D -->|No| F["Use zero endpoint second derivatives"]
    E --> G["Build local cubic coefficients"]
    F --> G
    G --> H{"Requested representation?"}
    H -->|Lagrange| I{"Query exactly a knot?"}
    I -->|Yes| J["Return recorded value"]
    I -->|No| K["Evaluate scaled barycentric ratio"]
    H -->|Newton| L["Evaluate nested polynomial"]
    H -->|Spline| M{"Query in knot range?"}
    M -->|No| N["Raise query-domain error"]
    M -->|Yes| O["Binary search interval; evaluate local cubic"]

The diagram assumes successful finite arithmetic and an accepted pivot sequence. An invalid input, zero computed barycentric denominator, unaccepted pivot, or nonfinite checked result raises an exception. Building the combined model prepares all three representations; an application needing only one can separate the preprocessing.

Worked Example with a Full Trace

Use $(x_i,y_i)=(0,0),(1,1),(2,0),(3,1)$ and query $q=1/2$. First compute every Lagrange basis contribution:

The basis values sum to $1$, and the contributions sum to $P(1/2)=1$. Raw barycentric weights are $[-1/6,1/2,-1/2,1/6]$; the normalized weights used in code are $[-1/3,1,-1,1/3]$. After scaling each term by the nearest-node distance $1/2$, the four denominator terms are $[-1/3,-1,1/3,-1/15]$. Both the numerator and denominator total $-16/15$, giving the same value $1$.

The complete divided-difference triangle is:

Thus $c=[0,1,-1,2/3]$ and $P(q)=q-q(q-1)+(2/3)q(q-1)(q-2)$. Nested evaluation at $1/2$ visits every coefficient:

1. Start $v=2/3$.
2. At index $2$, $v=-1+(1/2-2)(2/3)=-2$.
3. At index $1$, $v=1+(1/2-1)(-2)=2$.
4. At index $0$, $v=0+(1/2-0)(2)=1$.

For the natural spline, all interval widths are $1$ and the secant slopes are $[1,-1,1]$. With $M_0=M_3=0$, the complete interior system is

$$ \begin{bmatrix}4&1\\1&4\end{bmatrix} \begin{bmatrix}M_1\\M_2\end{bmatrix} =\begin{bmatrix}-12\\12\end{bmatrix}. $$

The first pivot is $4$; the maximum-magnitude row is already first. Eliminate the lower entry with multiplier $1/4$, obtaining the second row $[0,15/4]$ with right-hand side $15$. Back-substitution gives $M_2=4$, then $M_1=(-12-4)/4=-4$. The code first divides the whole system by $4$, so its pivot values are $1$ and $15/16$ with the same multiplier and final solution.

All three local spline pieces are now determined:

The query belongs to the first interval, so $S(1/2)=(5/3)(1/2)-(2/3)(1/2)^3=3/4$. At the interior knot $1$, the adjoining first derivatives both equal $-1/3$ and the second derivatives both equal $-4$; at knot $2$, they are $-1/3$ and $4$. The endpoint second derivatives are zero, completing the natural boundary conditions.

Implementation in Python, C++, and Java

Each complete program constructs a combined Models object from nested $[x,y]$ rows. It exposes lagrange(q), newton(q), and spline(q), along with normalized barycentric weights, Newton coefficients, spline second derivatives, and local cubic coefficients. Inputs are copied, and construction does not reorder or modify the caller's data.

The generic spline-system solver uses partial row pivoting and a normalized pivot tolerance of $10^{-12}$. The class uses direct $O(n^2)$ weight construction and checks for lost or nonfinite weights; it is intended for moderate knot counts and numerical scales. Common weight normalization and query-distance scaling do not cure every extreme spacing or conditioning problem.

The driver uses exactly the four knots and query from the full trace. It checks known between-knot values, the maximum knot residual across all representations, agreement of the two polynomial forms, and the interior spline-system residual. The local coefficient rows have order $[a_i,b_i,c_i,d_i]$ in powers of $q-x_i$.

import math


def finite(v):
    if not math.isfinite(v):
        raise ValueError("nonfinite input or arithmetic result")
    return v


def solve_pivoted(matrix, rhs, tol):
    n = len(rhs)
    if n == 0 or len(matrix) != n or any(len(row) != n for row in matrix):
        raise ValueError("nonempty square system required")
    if not math.isfinite(tol) or not 0.0 < tol < 1.0:
        raise ValueError("relative pivot tolerance must be in (0, 1)")
    scale = max(abs(finite(v)) for row in matrix for v in row)
    if scale == 0.0:
        raise ValueError("zero matrix")
    a = [[finite(v/scale) for v in row] for row in matrix]
    b = [finite(finite(v)/scale) for v in rhs]
    for k in range(n):
        pivot = max(range(k,n), key=lambda i: abs(a[i][k]))
        if abs(a[pivot][k]) <= tol:
            raise ValueError("pivot below relative tolerance")
        a[k], a[pivot] = a[pivot], a[k]
        b[k], b[pivot] = b[pivot], b[k]
        for i in range(k+1,n):
            factor = a[i][k]/a[k][k]
            a[i][k] = 0.0
            for j in range(k+1,n):
                a[i][j] = finite(a[i][j]-factor*a[k][j])
            b[i] = finite(b[i]-factor*b[k])
    answer = [0.0]*n
    for i in range(n-1,-1,-1):
        total = b[i]
        for j in range(i+1,n):
            total = finite(total-a[i][j]*answer[j])
        answer[i] = finite(total/a[i][i])
    return answer


class Models:
    def __init__(self, data, tol=1e-12):
        if len(data) < 2 or any(len(row) != 2 for row in data):
            raise ValueError("at least two [x, y] rows required")
        if not math.isfinite(tol) or not 0.0 < tol < 1.0:
            raise ValueError("relative pivot tolerance must be in (0, 1)")
        self.x = [finite(row[0]) for row in data]
        self.y = [finite(row[1]) for row in data]
        n = len(data)
        if any(self.x[i] <= self.x[i-1] for i in range(1,n)):
            raise ValueError("nodes must be strictly increasing")
        finite(self.x[-1]-self.x[0])
        self.w = []
        for i in range(n):
            weight = 1.0
            for j in range(n):
                if i != j:
                    weight = finite(weight/(self.x[i]-self.x[j]))
            self.w.append(weight)
        scale = max(abs(v) for v in self.w)
        if scale == 0.0:
            raise ValueError("barycentric weights underflowed")
        self.w = [v/scale for v in self.w]
        if any(v == 0.0 for v in self.w):
            raise ValueError("barycentric weight underflow")
        self.c = self.y.copy()
        for order in range(1,n):
            for i in range(n-1,order-1,-1):
                self.c[i] = finite((self.c[i]-self.c[i-1])/(self.x[i]-self.x[i-order]))
        h = [self.x[i+1]-self.x[i] for i in range(n-1)]
        slope = [finite((self.y[i+1]-self.y[i])/h[i]) for i in range(n-1)]
        self.second = [0.0]*n
        if n > 2:
            size = n-2
            a, rhs = [[0.0]*size for _ in range(size)], [0.0]*size
            for j in range(size):
                i = j+1
                a[j][j] = finite(2.0*(h[i-1]+h[i]))
                if j > 0:
                    a[j][j-1] = h[i-1]
                if j+1 < size:
                    a[j][j+1] = h[i]
                rhs[j] = finite(6.0*(slope[i]-slope[i-1]))
            self.second[1:-1] = solve_pivoted(a,rhs,tol)
        self.pieces = []
        for i in range(n-1):
            m0, m1 = self.second[i], self.second[i+1]
            b = finite(slope[i]-h[i]*(2.0*m0+m1)/6.0)
            self.pieces.append([self.y[i],b,m0/2.0,finite((m1-m0)/(6.0*h[i]))])

    def lagrange(self, q):
        finite(q)
        for x,y in zip(self.x,self.y):
            if q == x:
                return y
        distances = [finite(q-x) for x in self.x]
        scale = min(abs(d) for d in distances)
        numerator, denominator = 0.0, 0.0
        for w,d,y in zip(self.w,distances,self.y):
            term = w*(scale/d)
            numerator = finite(numerator+term*y)
            denominator = finite(denominator+term)
        if denominator == 0.0:
            raise ValueError("barycentric denominator canceled to zero")
        return finite(numerator/denominator)

    def newton(self, q):
        finite(q)
        value = self.c[-1]
        for i in range(len(self.c)-2,-1,-1):
            value = finite(self.c[i]+finite(q-self.x[i])*value)
        return value

    def spline(self, q):
        finite(q)
        if q < self.x[0] or q > self.x[-1]:
            raise ValueError("spline query outside knot interval")
        lo, hi = 0, len(self.x)-1
        while lo+1 < hi:
            mid = (lo+hi)//2
            if q < self.x[mid]:
                hi = mid
            else:
                lo = mid
        a,b,c,d = self.pieces[lo]
        t = q-self.x[lo]
        return finite(((d*t+c)*t+b)*t+a)


def main():
    data = [[0.0,0.0],[1.0,1.0],[2.0,0.0],[3.0,1.0]]
    model, q, tol = Models(data), 0.5, 1e-12
    p, n, s = model.lagrange(q), model.newton(q), model.spline(q)
    residual = max(abs(method(x)-y) for method in [model.lagrange,model.newton,model.spline]
                   for x,y in data)
    system = 0.0
    for i in range(1,len(data)-1):
        h0, h1 = data[i][0]-data[i-1][0], data[i+1][0]-data[i][0]
        rhs = 6.0*((data[i+1][1]-data[i][1])/h1-(data[i][1]-data[i-1][1])/h0)
        lhs = h0*model.second[i-1]+2*(h0+h1)*model.second[i]+h1*model.second[i+1]
        system = max(system,abs(lhs-rhs))
    ok = max(abs(p-1.0),abs(n-1.0),abs(s-0.75),residual,system) <= tol
    print(f"lagrange P(0.5) = {p:.6f}")
    print(f"newton P(0.5) = {n:.6f}")
    print(f"natural spline S(0.5) = {s:.6f}")
    print(f"max knot residual = {residual:.6f}")
    print(f"polynomial agreement error = {abs(p-n):.6f}")
    print(f"spline system residual = {system:.6f}")
    print(f"check = {'PASS' if ok else 'FAIL'}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>
using Matrix=std::vector<std::vector<double>>;
double finite_value(double v) {
    if(!std::isfinite(v)) throw std::invalid_argument("nonfinite input or arithmetic result");
    return v;
}
std::vector<double> solve_pivoted(const Matrix& matrix,const std::vector<double>& rhs,double tol) {
    int n=static_cast<int>(rhs.size());
    if(n==0 || static_cast<int>(matrix.size())!=n) throw std::invalid_argument("nonempty square system required");
    if(!std::isfinite(tol) || tol<=0 || tol>=1) throw std::invalid_argument("invalid relative pivot tolerance");
    double scale=0;
    for(const auto& row:matrix) {
        if(static_cast<int>(row.size())!=n) throw std::invalid_argument("nonempty square system required");
        for(double v:row) scale=std::max(scale,std::abs(finite_value(v)));
    }
    if(scale==0) throw std::invalid_argument("zero matrix");
    Matrix a=matrix;auto b=rhs;
    for(int i=0;i<n;++i) {
        for(double& v:a[i]) v=finite_value(v/scale);
        b[i]=finite_value(finite_value(b[i])/scale);
    }
    for(int k=0;k<n;++k) {
        int pivot=k;for(int i=k+1;i<n;++i) if(std::abs(a[i][k])>std::abs(a[pivot][k])) pivot=i;
        if(std::abs(a[pivot][k])<=tol) throw std::invalid_argument("pivot below relative tolerance");
        std::swap(a[k],a[pivot]);std::swap(b[k],b[pivot]);
        for(int i=k+1;i<n;++i) {
            double factor=a[i][k]/a[k][k];a[i][k]=0;
            for(int j=k+1;j<n;++j) a[i][j]=finite_value(a[i][j]-factor*a[k][j]);
            b[i]=finite_value(b[i]-factor*b[k]);
        }
    }
    std::vector<double> answer(n);
    for(int i=n-1;i>=0;--i) {
        double total=b[i];for(int j=i+1;j<n;++j) total=finite_value(total-a[i][j]*answer[j]);
        answer[i]=finite_value(total/a[i][i]);
    }
    return answer;
}
struct Models {
    std::vector<double> x,y,w,c,second;Matrix pieces;
    explicit Models(const Matrix& data,double tol=1e-12) {
        int n=static_cast<int>(data.size());
        if(n<2) throw std::invalid_argument("at least two [x, y] rows required");
        if(!std::isfinite(tol) || tol<=0 || tol>=1) throw std::invalid_argument("invalid relative pivot tolerance");
        for(const auto& row:data) {
            if(row.size()!=2) throw std::invalid_argument("[x, y] rows required");
            x.push_back(finite_value(row[0]));y.push_back(finite_value(row[1]));
        }
        for(int i=1;i<n;++i) if(x[i]<=x[i-1]) throw std::invalid_argument("nodes must be strictly increasing");
        finite_value(x.back()-x.front());w.resize(n);
        for(int i=0;i<n;++i) {
            double weight=1;for(int j=0;j<n;++j) if(i!=j) weight=finite_value(weight/(x[i]-x[j]));
            w[i]=weight;
        }
        double scale=0;for(double v:w) scale=std::max(scale,std::abs(v));
        if(scale==0) throw std::invalid_argument("barycentric weights underflowed");
        for(double& v:w) {v/=scale;if(v==0) throw std::invalid_argument("barycentric weight underflow");}
        c=y;
        for(int order=1;order<n;++order)
            for(int i=n-1;i>=order;--i) c[i]=finite_value((c[i]-c[i-1])/(x[i]-x[i-order]));
        std::vector<double> h(n-1),slope(n-1);second.assign(n,0);
        for(int i=0;i<n-1;++i) {h[i]=x[i+1]-x[i];slope[i]=finite_value((y[i+1]-y[i])/h[i]);}
        if(n>2) {
            int size=n-2;Matrix a(size,std::vector<double>(size));std::vector<double> rhs(size);
            for(int j=0;j<size;++j) {
                int i=j+1;a[j][j]=finite_value(2.0*(h[i-1]+h[i]));
                if(j>0) a[j][j-1]=h[i-1];
                if(j+1<size) a[j][j+1]=h[i];
                rhs[j]=finite_value(6.0*(slope[i]-slope[i-1]));
            }
            auto interior=solve_pivoted(a,rhs,tol);for(int i=0;i<size;++i) second[i+1]=interior[i];
        }
        for(int i=0;i<n-1;++i) {
            double m0=second[i],m1=second[i+1],b=finite_value(slope[i]-h[i]*(2.0*m0+m1)/6.0);
            pieces.push_back({y[i],b,m0/2.0,finite_value((m1-m0)/(6.0*h[i]))});
        }
    }
    double lagrange(double q) const {
        finite_value(q);int n=static_cast<int>(x.size());
        for(int i=0;i<n;++i) if(q==x[i]) return y[i];
        std::vector<double> distances(n);double scale=std::abs(finite_value(q-x[0]));
        for(int i=0;i<n;++i) {distances[i]=finite_value(q-x[i]);scale=std::min(scale,std::abs(distances[i]));}
        double numerator=0,denominator=0;
        for(int i=0;i<n;++i) {
            double term=w[i]*(scale/distances[i]);
            numerator=finite_value(numerator+term*y[i]);denominator=finite_value(denominator+term);
        }
        if(denominator==0) throw std::invalid_argument("barycentric denominator canceled to zero");
        return finite_value(numerator/denominator);
    }
    double newton(double q) const {
        finite_value(q);double value=c.back();
        for(int i=static_cast<int>(c.size())-2;i>=0;--i) value=finite_value(c[i]+finite_value(q-x[i])*value);
        return value;
    }
    double spline(double q) const {
        finite_value(q);if(q<x.front() || q>x.back()) throw std::invalid_argument("spline query outside knot interval");
        int lo=0,hi=static_cast<int>(x.size())-1;
        while(lo+1<hi) {int mid=(lo+hi)/2;if(q<x[mid]) hi=mid;else lo=mid;}
        const auto& p=pieces[lo];double t=q-x[lo];return finite_value(((p[3]*t+p[2])*t+p[1])*t+p[0]);
    }
};
int main() {
    Matrix data={ {0.,0.}, {1.,1.}, {2.,0.}, {3.,1.} };Models model(data);
    double q=.5,tol=1e-12,p=model.lagrange(q),n=model.newton(q),s=model.spline(q),residual=0,system=0;
    for(const auto& row:data) {
        residual=std::max(residual,std::abs(model.lagrange(row[0])-row[1]));
        residual=std::max(residual,std::abs(model.newton(row[0])-row[1]));
        residual=std::max(residual,std::abs(model.spline(row[0])-row[1]));
    }
    for(int i=1;i<static_cast<int>(data.size())-1;++i) {
        double h0=data[i][0]-data[i-1][0],h1=data[i+1][0]-data[i][0];
        double rhs=6.0*((data[i+1][1]-data[i][1])/h1-(data[i][1]-data[i-1][1])/h0);
        double lhs=h0*model.second[i-1]+2*(h0+h1)*model.second[i]+h1*model.second[i+1];
        system=std::max(system,std::abs(lhs-rhs));
    }
    bool ok=std::max({std::abs(p-1.),std::abs(n-1.),std::abs(s-.75),residual,system})<=tol;
    std::cout<<std::fixed<<std::setprecision(6);
    std::cout<<"lagrange P(0.5) = "<<p<<'\n'<<"newton P(0.5) = "<<n<<'\n';
    std::cout<<"natural spline S(0.5) = "<<s<<'\n'<<"max knot residual = "<<residual<<'\n';
    std::cout<<"polynomial agreement error = "<<std::abs(p-n)<<'\n'<<"spline system residual = "<<system<<'\n';
    std::cout<<"check = "<<(ok?"PASS":"FAIL")<<'\n';return 0;
}

import java.util.Locale;
public class Main {
    static double finite(double v) {
        if(!Double.isFinite(v)) throw new IllegalArgumentException("nonfinite input or arithmetic result");
        return v;
    }
    static double[] solvePivoted(double[][] matrix,double[] rhs,double tol) {
        int n=rhs.length;
        if(n==0 || matrix.length!=n) throw new IllegalArgumentException("nonempty square system required");
        if(!Double.isFinite(tol) || tol<=0 || tol>=1) throw new IllegalArgumentException("invalid relative pivot tolerance");
        double scale=0;
        for(double[] row:matrix) {
            if(row.length!=n) throw new IllegalArgumentException("nonempty square system required");
            for(double v:row) scale=Math.max(scale,Math.abs(finite(v)));
        }
        if(scale==0) throw new IllegalArgumentException("zero matrix");
        double[][] a=new double[n][n];double[] b=new double[n];
        for(int i=0;i<n;++i) {for(int j=0;j<n;++j) a[i][j]=finite(matrix[i][j]/scale);b[i]=finite(finite(rhs[i])/scale);}
        for(int k=0;k<n;++k) {
            int pivot=k;for(int i=k+1;i<n;++i) if(Math.abs(a[i][k])>Math.abs(a[pivot][k])) pivot=i;
            if(Math.abs(a[pivot][k])<=tol) throw new IllegalArgumentException("pivot below relative tolerance");
            double[] row=a[k];a[k]=a[pivot];a[pivot]=row;double t=b[k];b[k]=b[pivot];b[pivot]=t;
            for(int i=k+1;i<n;++i) {
                double factor=a[i][k]/a[k][k];a[i][k]=0;
                for(int j=k+1;j<n;++j) a[i][j]=finite(a[i][j]-factor*a[k][j]);
                b[i]=finite(b[i]-factor*b[k]);
            }
        }
        double[] answer=new double[n];
        for(int i=n-1;i>=0;--i) {
            double total=b[i];for(int j=i+1;j<n;++j) total=finite(total-a[i][j]*answer[j]);
            answer[i]=finite(total/a[i][i]);
        }
        return answer;
    }
    static class Models {
        double[] x,y,w,c,second;double[][] pieces;
        Models(double[][] data,double tol) {
            int n=data.length;
            if(n<2) throw new IllegalArgumentException("at least two [x, y] rows required");
            if(!Double.isFinite(tol) || tol<=0 || tol>=1) throw new IllegalArgumentException("invalid relative pivot tolerance");
            x=new double[n];y=new double[n];w=new double[n];
            for(int i=0;i<n;++i) {
                if(data[i].length!=2) throw new IllegalArgumentException("[x, y] rows required");
                x[i]=finite(data[i][0]);y[i]=finite(data[i][1]);
            }
            for(int i=1;i<n;++i) if(x[i]<=x[i-1]) throw new IllegalArgumentException("nodes must be strictly increasing");
            finite(x[n-1]-x[0]);
            for(int i=0;i<n;++i) {
                double weight=1;for(int j=0;j<n;++j) if(i!=j) weight=finite(weight/(x[i]-x[j]));w[i]=weight;
            }
            double scale=0;for(double v:w) scale=Math.max(scale,Math.abs(v));
            if(scale==0) throw new IllegalArgumentException("barycentric weights underflowed");
            for(int i=0;i<n;++i) {w[i]/=scale;if(w[i]==0) throw new IllegalArgumentException("barycentric weight underflow");}
            c=y.clone();
            for(int order=1;order<n;++order)
                for(int i=n-1;i>=order;--i) c[i]=finite((c[i]-c[i-1])/(x[i]-x[i-order]));
            double[] h=new double[n-1],slope=new double[n-1];second=new double[n];
            for(int i=0;i<n-1;++i) {h[i]=x[i+1]-x[i];slope[i]=finite((y[i+1]-y[i])/h[i]);}
            if(n>2) {
                int size=n-2;double[][] a=new double[size][size];double[] rhs=new double[size];
                for(int j=0;j<size;++j) {
                    int i=j+1;a[j][j]=finite(2.0*(h[i-1]+h[i]));
                    if(j>0) a[j][j-1]=h[i-1];
                    if(j+1<size) a[j][j+1]=h[i];
                    rhs[j]=finite(6.0*(slope[i]-slope[i-1]));
                }
                double[] interior=solvePivoted(a,rhs,tol);System.arraycopy(interior,0,second,1,size);
            }
            pieces=new double[n-1][4];
            for(int i=0;i<n-1;++i) {
                double m0=second[i],m1=second[i+1],b=finite(slope[i]-h[i]*(2.0*m0+m1)/6.0);
                pieces[i]=new double[]{y[i],b,m0/2.0,finite((m1-m0)/(6.0*h[i]))};
            }
        }
        double lagrange(double q) {
            finite(q);int n=x.length;for(int i=0;i<n;++i) if(q==x[i]) return y[i];
            double[] distances=new double[n];double scale=Math.abs(finite(q-x[0]));
            for(int i=0;i<n;++i) {distances[i]=finite(q-x[i]);scale=Math.min(scale,Math.abs(distances[i]));}
            double numerator=0,denominator=0;
            for(int i=0;i<n;++i) {
                double term=w[i]*(scale/distances[i]);numerator=finite(numerator+term*y[i]);denominator=finite(denominator+term);
            }
            if(denominator==0) throw new IllegalArgumentException("barycentric denominator canceled to zero");
            return finite(numerator/denominator);
        }
        double newton(double q) {
            finite(q);double value=c[c.length-1];
            for(int i=c.length-2;i>=0;--i) value=finite(c[i]+finite(q-x[i])*value);
            return value;
        }
        double spline(double q) {
            finite(q);if(q<x[0] || q>x[x.length-1]) throw new IllegalArgumentException("spline query outside knot interval");
            int lo=0,hi=x.length-1;while(lo+1<hi) {int mid=(lo+hi)/2;if(q<x[mid]) hi=mid;else lo=mid;}
            double[] p=pieces[lo];double t=q-x[lo];return finite(((p[3]*t+p[2])*t+p[1])*t+p[0]);
        }
    }
    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);double[][] data={ {0.,0.}, {1.,1.}, {2.,0.}, {3.,1.} };
        Models model=new Models(data,1e-12);double q=.5,tol=1e-12,p=model.lagrange(q),n=model.newton(q),s=model.spline(q),residual=0,system=0;
        for(double[] row:data) {
            residual=Math.max(residual,Math.abs(model.lagrange(row[0])-row[1]));
            residual=Math.max(residual,Math.abs(model.newton(row[0])-row[1]));
            residual=Math.max(residual,Math.abs(model.spline(row[0])-row[1]));
        }
        for(int i=1;i<data.length-1;++i) {
            double h0=data[i][0]-data[i-1][0],h1=data[i+1][0]-data[i][0];
            double rhs=6.0*((data[i+1][1]-data[i][1])/h1-(data[i][1]-data[i-1][1])/h0);
            double lhs=h0*model.second[i-1]+2*(h0+h1)*model.second[i]+h1*model.second[i+1];system=Math.max(system,Math.abs(lhs-rhs));
        }
        boolean ok=Math.max(Math.max(Math.abs(p-1.),Math.abs(n-1.)),Math.max(Math.abs(s-.75),Math.max(residual,system)))<=tol;
        System.out.printf("lagrange P(0.5) = %.6f%nnewton P(0.5) = %.6f%n",p,n);
        System.out.printf("natural spline S(0.5) = %.6f%nmax knot residual = %.6f%n",s,residual);
        System.out.printf("polynomial agreement error = %.6f%nspline system residual = %.6f%n",Math.abs(p-n),system);
        System.out.println("check = "+(ok?"PASS":"FAIL"));
    }
}

All three print:

• lagrange P(0.5) = 1.000000
• newton P(0.5) = 1.000000
• natural spline S(0.5) = 0.750000
• max knot residual = 0.000000
• polynomial agreement error = 0.000000
• spline system residual = 0.000000
• check = PASS

A knot residual can be zero by construction; barycentric evaluation deliberately returns the supplied value at an exact knot. That check alone cannot validate all coefficients or establish between-knot accuracy. The known query values and spline system give additional checks; derivative-continuity and independent-reference tests are useful when changing the implementation.

Complexity, Correctness, and Invariants

Let $n$ be the number of knots. Separate preprocessing from query cost: an interpolant may be built once and evaluated thousands of times. The space column distinguishes stored model data from temporary construction workspace.

The combined class builds all three models, so its dense solve dominates preprocessing and peak space for large $n$. With two knots no system is solved. For ordinary valid inputs of a fixed size, the construction loop bounds are data-independent; pivot choices affect arithmetic and row swaps but not the stated asymptotic order. Newton query work is fixed, while barycentric exact-knot exits and binary-search paths vary. A distribution-free average query cost is not defined without a query model. If a spline's interval index is already known, its evaluation itself is $O(1)$.

Polynomial existence and uniqueness. Lagrange's cardinal identities give a polynomial meeting every knot condition. If two polynomials of degree at most $n-1$ met those conditions, their difference would have $n$ distinct roots while having degree at most $n-1$, so it must be zero. Therefore the Lagrange and Newton forms represent the same exact polynomial.

Newton's construction invariant. At divided-difference order $k$, the entries about to be used belong to order $k-1$. Descending index updates preserve those inputs until they have been consumed; the diagonal entry at index $k$ becomes the coefficient for $\prod_{j=0}^{k-1}(q-x_j)$. Nested evaluation maintains the remaining suffix of that polynomial without expanding its powers.

Spline existence and uniqueness. Each local cubic is determined by its endpoint values and second derivatives. The interior equations impose continuity of the first derivative; shared $M_i$ impose continuity of the second derivative. With natural endpoint values fixed, the interior matrix has positive diagonal, positive interval widths, and strict diagonal dominance. Symmetry and positivity give a unique solution. Partial pivoting is an equivalent row-operation solve; the normalized tolerance can still reject a numerically difficult representation even when the mathematical system is nonsingular.

If the data are exact samples of a sufficiently smooth function $f$, the polynomial interpolation error satisfies

$$ f(q)-P(q)=\frac{f^{(n)}(\xi)}{n!}\prod_{i=0}^{n-1}(q-x_i) $$

for some $\xi$ in the smallest interval containing the nodes and query. This is not an available numerical bound unless the derivative is bounded independently. The nodal product explains why node placement matters. Spline convergence bounds depend on mesh regularity, smoothness, and compatible boundary information; a natural endpoint condition should not be assumed to give a fourth-order global approximation for every smooth function with nonzero endpoint curvature.

Variants, Pitfalls, and Optimizations

Avoid confusing representation with the interpolation model. Switching from direct Lagrange products to barycentric Lagrange can improve evaluation, but it leaves the polynomial unchanged. Switching to a natural spline changes the model. A stable evaluator cannot eliminate oscillation inherent in a chosen high-degree interpolating polynomial.

Runge's example makes this visible. For $f(x)=1/(1+25x^2)$, use $11$ equally spaced nodes on $[-1,1]$. At $q=0.94$, the degree-at-most-ten polynomial is about $1.958952$, while the function value is about $0.043309$. The natural spline is about $0.043824$. All interpolants still meet their knot values. The large polynomial error near the endpoints is a node-placement and model issue, not disagreement between Lagrange and Newton in exact arithmetic.

If node placement is controllable, Chebyshev-distributed nodes can reduce the endpoint oscillation associated with equally spaced high-degree interpolation. Piecewise low-degree methods are another option. If nodes are fixed measurements, changing them is not simply an implementation detail.

Specify boundary conditions. Natural means zero endpoint second derivative, not zero first derivative or a guarantee of a straight endpoint interval. Clamped splines specify endpoint slopes; periodic splines join matching endpoint behavior; not-a-knot conditions select another spline. In particular, a natural spline does not reproduce every quadratic or cubic: the natural interpolant of $x^2$ at $0,1,2$ gives $S(1/2)=5/16$, rather than $1/4$, because the quadratic's endpoint curvature is not zero. When comparing with a library, explicitly select natural boundary conditions; for example, SciPy's CubicSpline default is not-a-knot.

Do not assume shape preservation. Ordinary cubic splines can overshoot monotone data or become negative between positive knots. Piecewise linear interpolation preserves interval endpoint bounds, while a shape-preserving cubic method such as PCHIP may be a better choice when monotonicity matters. A smoothing spline or regression is often more appropriate when data contain noise: forcing exact interpolation can reproduce noise as if it were signal.

Handle duplicates and scaling deliberately. Equal nodes with unequal values cannot satisfy a single-valued interpolation function. Equal nodes with derivative information require a Hermite formulation; they are not valid denominators for these routines. Centering and scaling the horizontal coordinate can reduce arithmetic problems. Clustered nodes may also make the interpolation problem sensitive even after scaling. A floating-point pivot threshold belongs to the solve policy; it should not silently redefine distinct knots as identical.

Improve repeated work. Retain barycentric weights when only the $y$ values change, reuse the spline-system factorization when the same knots support multiple data channels, and exploit the tridiagonal system for linear-time spline construction. Newton's form supports adding a distinct knot with an additional coefficient in $O(n)$ arithmetic if an incremental interface is maintained; the supplied combined class rebuilds its representations instead. For sorted spline queries, advance through intervals once rather than binary-searching every query independently.

Treat extrapolation as a separate decision. A polynomial can be evaluated beyond its nodes, but its magnitude and sensitivity may grow rapidly. Natural boundary conditions do not justify linear extrapolation outside the domain. The spline API here rejects such queries; a different application must deliberately choose and document an extension rule.

Practice Problems and Real-World Use

Begin with exact polynomials and small rational data, then test nonuniform spacing, boundary conditions, and adverse examples. In addition to knot residuals, compare the two polynomial forms and check both sides of every spline join for value, slope, and curvature continuity.

For further reading, NIST's interpolation chapter develops Lagrange, barycentric, and Newton formulas and their remainder terms. SciPy's cubic spline documentation states the boundary-condition choices, while its barycentric interpolator documentation discusses polynomial representation and node-placement concerns. Match the mathematical model and boundary conditions before comparing numerical outputs.