Natural-cubic tridiagonal Cholesky solve with a residual certificate

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 ↗.