Spline evaluation and interpolation algorithms

Use zero-based indices for this procedure: write the clamped knot vector from B-spline spaces and knot conventions as $t_0\le \cdots \le t_{n+p}$, where $t_i={\xi }_{i+1}$, with $n\ge p+1$, and write the real degree-$p$ coefficients as $c_0,\dots ,c_{n-1}$. The input includes $p\in {\mathbb{N}}_0$, a nonempty span index $k$ with $p\le k\le n-1$ and $t_k<t_{k+1}$, a real query, and the $p+1$ active coefficients $c_{k-p},\dots ,c_k$. The query is either $x\in [t_k,t_{k+1})$ or, for the left limit at the right endpoint $b$, $x=b$ with the last nonempty span supplied. At an interior knot, the half-open convention selects the span on its right. In exact mode the knots, query, and coefficients are real numbers. In floating-point mode the knots, query, and active coefficients are finite stored real numbers in the selected format; conversion into that format is outside this procedure.

The arithmetic mode is either exact real arithmetic or the operationwise floating-point model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. The output is success(value), invalid-input, or, in floating-point mode, range-failure(operation, r, j). The span search is supplied by the caller and is not part of the evaluation loop. Invalid indices, coefficient counts, span or query conditions, or floating-point input formats return invalid-input.

1. If the indices, coefficient count, span, or query do not satisfy the input conditions, return invalid-input.
2. If $p=0$, return $c_k$ by copy.
3. Initialize $d_j^{(0)}=c_{k-p+j}$ for $0\le j\le p$.
4. For $r=1,\dots ,p$, and then $j=p,p-1,\dots ,r$, put $L=k-p+j$ and $U=k+1+j-r$. In exact mode set $\alpha =(x-t_L)\mathrm{/}(t_U-t_L)$ and $d_j^{(r)}=(1-\alpha )d_{j-1}^{(r-1)}+\alpha d_j^{(r-1)}$. In floating-point mode first compute $v=\mathrm{fl}(x-t_L)$. On failure return range-failure for the numerator subtraction. If $v=0$, copy $d_{j-1}^{(r-1)}$ to $d_j^{(r)}$ and skip the denominator subtraction, division, complement subtraction, products, and addition. Otherwise compute $w=\mathrm{fl}(t_U-t_L)$ and then $\widehat{\alpha }=\mathrm{fl}(v\mathrm{/}w)$; on the first failure return its operation and indices. Next compute $\widehat{\beta }=\mathrm{fl}(1-\widehat{\alpha })$, followed by $a=\mathrm{fl}(\widehat{\beta }{d}_{j-1}^{(r-1)})$, $b=\mathrm{fl}(\widehat{\alpha }{d}_j^{(r-1)})$, and $d_j^{(r)}=\mathrm{fl}(a+b)$, in that order, returning the first range-failure. Keep the old $d_j^{(r-1)}$ until the update is complete.
5. Return success(d_p^{(p)}).

For every loop entry, $L\le k<U$, so $t_L\le t_k<t_{k+1}\le t_U$ and the denominator is positive. The exact weight therefore lies in $[0,1]$. When a span endpoint is requested as a one-sided value, the same recurrence is evaluated at that endpoint with the specified half-open or left-limit convention.

Use the zero-based reindexing $t_i={\xi }_{i+1}$ of the clamped knot vector in B-spline spaces and knot conventions, so $t_0\le \cdots \le t_{n+p}$ with $n\ge p+1$. Let $p\in {\mathbb{N}}_0$, $d\in \mathbb{N}$, and let the controls be $P_0,\dots ,P_{n-1}\in {\mathbb{R}}^d$. Let $k$ be a nonempty span index with $t_k<u<t_{k+1}$, where $u$ is not already a knot. Insert $u$ between $t_k$ and $t_{k+1}$ to obtain $t^{\mathrm{\prime }}$. The output controls $Q_0,\dots ,Q_n$ are defined by

$Q_i=\{\begin{array}{ll}{P}_i,& 0\le i\le k-p,\\ {\alpha }_i{P}_i+(1-{\alpha }_i)P_{i-1},& k-p+1\le i\le k,\\ P_{i-1},& k+1\le i\le n,\end{array}{\alpha }_i=\frac{u-t_i}{t_{i+p}-t_i}\mathrm{.}$

The middle range is empty when $p=0$. The curve uses the same degree and the B-spline basis for the new knot vector. Each denominator in the middle range is positive and each ${\alpha }_i$ lies in $[0,1]$, so the update is an exact affine blend. This procedure inserts one knot that is new to the vector; repeated insertion at an existing knot is not part of this algorithm's input contract.

Let $a=x_0<x_1<\cdots <x_m=b$, with $m\ge 1$, let $y_0,\dots ,y_m\in \mathbb{R}$, and prescribe endpoint slopes ${\sigma }_a,{\sigma }_b\in \mathbb{R}$. A clamped cubic spline interpolant for these data is a function $s\in C^2([a,b])$ whose restriction to every $[x_i,x_{i+1}]$ is a polynomial of degree at most three, which satisfies $s(x_i)=y_i$ for every $i$, and which satisfies $s^{\mathrm{\prime }}(a)={\sigma }_a$ and $s^{\mathrm{\prime }}(b)={\sigma }_b$.

Here $C^2$ means that the function and its first two derivatives are continuous on $[a,b]$, with one-sided derivatives at the endpoints. This uses the piecewise-cubic convention in Natural cubic splines and replaces its natural endpoint-curvature conditions by prescribed endpoint slopes. It does not alter the natural cubic interpolant or its system.

Let $x_1<\cdots <x_r$ be natural interpolation knots, $r\ge 2$, and put $h_i=x_{i+1}-x_i$ and $H={\max }_i{h}_i$. Let $s$ be the natural cubic interpolant with data $y_i$ and nodal second derivatives $M_i$, with $M_1=M_r=0$. Let ${\widehat{M}}_1={\widehat{M}}_r=0$ and let $\widehat{s}$ be the piecewise cubic formed from the same values $y_i$ and the second derivatives ${\widehat{M}}_i$ using the natural cubic piece formula of Natural cubic interpolation at distinct knots. On the knot hull $[x_1,x_r]$, if $E={\max }_{1\le i\le r}\mathrm{\mid }{\widehat{M}}_i-M_i\mathrm{\mid },$ then $\mathrm{\parallel }\widehat{s}-s{\mathrm{\parallel }}_{L^{\mathrm{\infty }}([x_1,x_r])}\le \frac{H^2}{8}E\mathrm{.}$

Use one binary floating-point format from Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add with round-to-nearest, ties-to-even, a known rational unit roundoff $0<u\le 1\mathrm{/}4$, and the model's normal-or-zero success rule. Inputs are $r\ge 2$ finite stored real knots $x_1<\cdots <x_r$ and finite stored real values $y_1,\dots ,y_r$. These stored numbers are the exact interpolation data; no unrecorded perturbation of the data is included. Write $h_i=x_{i+1}-x_i$ for their exact rational differences. If a constant $2$ or $6$ needed by the selected branch is not representable as a finite normal number, return input-format-failure. Arithmetic statuses below report the first failing operation in the specified order. This algorithm makes no claim that floating-point success is automatic.

If $r=2$, return success with endpoint second derivatives ${\widehat{M}}_1={\widehat{M}}_2=0$ and error certificate $B=0$. The corresponding piece is the exact straight line. For $r\ge 3$, put $n=r-2$ and return certificate-parameter-failure unless $(n-1)u<1$.

For $r\ge 3$, assemble stored widths and secant slopes in increasing $i=1,\dots ,r-1$: compute ${\widehat{h}}_i=\mathrm{fl}(x_{i+1}-x_i)$, $a_i=\mathrm{fl}(y_{i+1}-y_i)$, and $q_i=\mathrm{fl}(a_i\mathrm{/}{\widehat{h}}_i)$, in that order. At the first failed operation return assembly-range-failure(i,operation). Successful width subtraction gives ${\widehat{h}}_i>0$. In increasing $k=1,\dots ,n$, for the row corresponding to knot $x_{k+1}$ compute

in that order. Copy $c_k={\widehat{h}}_{k+1}$ for $k<n$. On failure return assembly-range-failure(k,operation). Let $T\in {\mathbb{R}}^{n\times n}$ have diagonal $d_k$, adjacent off-diagonal entries $c_k$, and zero entries elsewhere; represent it by these three diagonals, not a dense array. Let the stored right-hand side be $f=(f_1,\dots ,f_n{)}^T$.

Factor $T$ with a real lower bidiagonal factor, stored as diagonal ${\mathrm{\ell }}_1,\dots ,{\mathrm{\ell }}_n$ and subdiagonal $e_1,\dots ,e_{n-1}$. For $k=1,\dots ,n$ in order, set $p_1=d_1$ by exact copy; for $k>1$ compute ${\zeta }_k=\mathrm{fl}(e_{k-1}{e}_{k-1})$ and $p_k=\mathrm{fl}(d_k-{\zeta }_k)$. If either operation fails, return factor-range-failure(k,operation). If the resulting $p_k$ fails the published positive-normal pivot test of Failure-aware rounded Cholesky factorization, return invalid-pivot(k,p_k). Otherwise compute ${\mathrm{\ell }}_k=\mathrm{fl}(\sqrt{p_k})$; on failure or a nonpositive or nonnormal result return factor-range-failure(k,operation). For $k<n$ compute $e_k=\mathrm{fl}(c_k\mathrm{/}{\mathrm{\ell }}_k)$ and return the same indexed factor status on failure.

Solve $Lz=f$ by forward substitution in increasing index order: set $z_1=\mathrm{fl}(f_1\mathrm{/}{\mathrm{\ell }}_1)$; for $k>1$, compute $a=\mathrm{fl}(e_{k-1}{z}_{k-1})$, $b=\mathrm{fl}(f_k-a)$, and $z_k=\mathrm{fl}(b\mathrm{/}{\mathrm{\ell }}_k)$. Solve $L^T\widehat{m}=z$ in decreasing order: set ${\widehat{m}}_n=\mathrm{fl}(z_n\mathrm{/}{\mathrm{\ell }}_n)$; for $k<n$, compute $a=\mathrm{fl}(e_k{\widehat{m}}_{k+1})$, $b=\mathrm{fl}(z_k-a)$, and ${\widehat{m}}_k=\mathrm{fl}(b\mathrm{/}{\mathrm{\ell }}_k)$. At the first failed operation return solve-range-failure(direction,k,operation). These are separate multiply, subtract, and divide operations; no fused multiply-add is used. Set ${\widehat{M}}_1={\widehat{M}}_r=0$ and copy the interior vector ${\widehat{M}}_{k+1}={\widehat{m}}_k$.

Compute the rounded residual in increasing row order. Within row $k$, form the products $T_{kj}{\widehat{m}}_j$ for existing columns $j$ in increasing column order, and add them sequentially from an exact zero accumulator to obtain ${\widehat{v}}_k$. Then compute ${\widehat{w}}_k=\mathrm{fl}({\widehat{v}}_k-f_k)$. Return the first residual-range-failure(k,j,operation) if any operation fails. Structural zero entries are omitted; all included products, additions, and the final subtraction use the operationwise model.

The certificate calculations that follow are exact rational arithmetic on the finite binary values and the rational $u$; they are not rounded floating-point calculations. Define

and

Form the exact natural-cubic matrix $A$ and right-hand side $b$ for the original stored data as in Natural cubic interpolation at distinct knots: $A_{kk}=2(h_k+h_{k+1})$, $A_{k,k+1}=A_{k+1,k}=h_{k+1}$, and $b_k=6((y_{k+2}-y_{k+1})\mathrm{/}{h}_{k+1}-(y_{k+1}-y_k)\mathrm{/}{h}_k)$, with only existing interior entries included; represent $A$ by its three diagonals and compute the following row quantities as a stream, without allocating dense matrices. Compute

and

All denominators are positive: $h_i>0$, $\delta >0$, and $1-(n-1)u>0$. Return success(\widehat M,C,\alpha,\beta,B). If a resource-limited implementation cannot complete the exact rational certificate, it may return certificate-resource-failure, which carries no certified $B$.

For $r\ge 3$, associate to the returned vector $\widehat{M}$ the piecewise cubic reconstruction on each $[x_i,x_{i+1}]$ given by

This reconstruction interpolates the stored data and need not be $C^1$ at an interior knot; it is a piecewise cubic reconstruction, not necessarily a natural spline. The certificate concerns exact evaluation of these pieces. Any floating-point evaluation of them needs a separate piece-evaluation error allowance.

This local procedure specializes the published rounded Cholesky and substitution methods to the natural-cubic tridiagonal system; it does not redefine natural interpolation. Its termination and linear arithmetic and storage bounds are proved in A successful natural-cubic tridiagonal Cholesky solve has a certified residual ↗.

Let $t_0\le \cdots \le t_{n+p}$ be a clamped degree-$p$ knot vector as in B-spline spaces and knot conventions, with endpoints $a=t_0$ and $b=t_{n+p}$. For a zero-based knot vector $v$, write $\xi (v{)}_j=v_{j-1}$ and $B_{i,q,v}=B_{i+1,q,\xi (v)}$. Let $k$ satisfy $t_k<t_{k+1}$, let $t_k<u<t_{k+1}$, and obtain $t^{\mathrm{\prime }}$ by inserting $u$ between $t_k$ and $t_{k+1}$. For $0\le q\le p$, define ${\lambda }_i^{(q)}=\{\begin{array}{ll}1,& i\le k-q,\\ \frac{u-t_i}{t_{i+q}-t_i},& k-q+1\le i\le k,\\ 0,& i\ge k+1.\end{array}$ Then for every degree $q\in \{0,\dots ,p\}$, every old basis index $0\le i\le n+p-q-1$, and every $x\in [a,b]$, $B_{i,q,t}(x)={\lambda }_i^{(q)}{B}_{i,q,t^{\mathrm{\prime }}}(x)+(1-{\lambda }_{i+1}^{(q)})B_{i+1,q,t^{\mathrm{\prime }}}(x)\mathrm{.}$

Let $p\in {\mathbb{N}}_0$, let $t_0\le \cdots \le t_{n+p}$ be the zero-based reindexing of a clamped knot vector, and let $S(x)={\sum }_{i=0}^{n-1}{c}_i{B}_{i,p}(x)$ be its degree-$p$ B-spline expansion, with $B_{i,p}=B_{i+1,p,\xi }$ under B-spline spaces and knot conventions. If $k$ is a nonempty span index and $x\in [t_k,t_{k+1})$, then the exact-arithmetic mode of de Boor coefficient evaluation on a specified knot span returns $S(x)$. For $x=b$ it returns the corresponding left limit when the last nonempty span is supplied.

For a supplied nonempty knot span and degree $p$, a successful call to the loop in de Boor coefficient evaluation on a specified knot span performs exactly $p(p+1)\mathrm{/}2$ coefficient updates when $p\ge 1$ and none when $p=0$. A call that returns a range failure performs no more than this many updates. Each update uses at most a constant number of scalar arithmetic operations, so the arithmetic work is $O(p^2)$ and the coefficient storage is $O(p)$; the caller's span search is excluded.

Let $a=x_0<x_1<\cdots <x_m=b$, $m\ge 1$, let $y_0,\dots ,y_m\in \mathbb{R}$, and let ${\sigma }_a,{\sigma }_b\in \mathbb{R}$. Put $h_i=x_{i+1}-x_i>0$ and $d_i=(y_{i+1}-y_i)\mathrm{/}{h}_i$. There is exactly one clamped cubic spline interpolant in the sense of Clamped cubic spline interpolant. Its nodal second derivatives $M_i=s^{\mathrm{\prime }\mathrm{\prime }}(x_i)$ are the unique solution of $\begin{array}{rl}{\displaystyle 2h_0{M}_0+h_0{M}_1}& {\displaystyle =6(d_0-{\sigma }_a),}\\ {\displaystyle h_{i-1}{M}_{i-1}+2(h_{i-1}+h_i)M_i+h_i{M}_{i+1}}& {\displaystyle =6(d_i-d_{i-1}),1\le i<m,}\\ {\displaystyle h_{m-1}{M}_{m-1}+2h_{m-1}{M}_m}& {\displaystyle =6({\sigma }_b-d_{m-1})\mathrm{.}}\end{array}$ On $[x_i,x_{i+1}]$ the interpolant is $s_i(x)=\frac{M_i(x_{i+1}-x{)}^3+M_{i+1}(x-x_i{)}^3}{6h_i}+(y_i-\frac{M_i{h}_i^2}{6})\frac{x_{i+1}-x}{h_i}+(y_{i+1}-\frac{M_{i+1}{h}_i^2}{6})\frac{x-x_i}{h_i}\mathrm{.}$