Finite-precision direct factorizations and least squares

At each precision, let the representable set, normal range, and unit roundoff be an instance of The standard relative floating-point model and unit roundoff specialized to a finite binary format of integer precision $p\ge 2$ and integer normal-exponent bounds $j_{\min }\le j\le j_{\max }$. Its normal values are $\pm m{2}^{j-p+1}$ for integers $2^{p-1}\le m<2^p$ and $j_{\min }\le j\le j_{\max }$, with $N_{\min }=2^{j_{\min }}$ and $N_{\max }=(2^p-1)2^{j_{\max }-p+1}$. The rounding rule is round-to-nearest, ties-to-even by the parity of the normalized significand. Each precision in a mixed-precision procedure has its own format, range, and valid unit-roundoff bound. On every real exact result that is zero or normal, $\mathrm{fl}$ uses this rule; every other nonzero exact result (subnormal or overflow), NaN, infinity, invalid operation, or division by zero returns range-failure. A real add, subtract, multiply, or divide is one rounded operation. For positive in-range $x$, $\mathrm{fl}(\sqrt{x})$ is correctly rounded to nearest, ties-to-even, and fails unless its result is positive and normal. A real fused multiply-add returns the correctly rounded value of the single exact expression $ab+c$; it is not two operations.

Complex addition and subtraction round real and imaginary components separately. A complex product (a+ib)(c+id) uses four rounded real products ac, bd, ad, bc, then rounded subtraction and addition, for six real operations. Conjugation and sign changes are exact. Scaling and division by a positive real scalar round each complex component separately; no other complex division primitive is used here.

A procedure succeeds only when each exact intermediate is zero or normal and every rounded output is finite; otherwise it returns the first range-failure with its operation or pivot index. For ku<1 write gamma_k=ku/(1-ku). This is a no-underflow success model, not the gradual-underflow model.

Fix an initialized finite state and an ordered list of $M\ge 0$ instructions, followed by a specified completion return. Each instruction reads only initialized data and deterministically either produces the next finite state or returns a specified early status in finite time. There are no backward jumps. Then execution terminates at the first early return, if any, or at the completion return, after at most $M$ attempted instructions. A finite nesting of counted loops with bounds fixed by the finite input has the same property whenever unrolling the loops gives such a list.

In particular, consider $d\ge 0$ levels containing $2^{d-\mathrm{\ell }}$ addition instructions at level $\mathrm{\ell }=1,\dots ,d$, with each addition returning a rounded value or immediate failure. A successful execution attempts exactly $2^d-1$ additions; every execution attempts at most $2^d-1$. For $d=0$ there are no additions. These conclusions assert finite execution, not arithmetic success or accuracy.

The source supplies background on rounded arithmetic and range limitations. The local argument proves the finite-schedule and counting assertions; it does not assume that a primitive operation can run indefinitely while searching for a certificate or a rounded answer.

Input real floating-point values $x_1,\dots ,x_n$, $n\ge 1$. Set $d=\lceil {\log }_{2}n\rceil$, place these values in leaves $1,\dots ,n$ from left to right, and append $2^d-n$ exact zero leaves. At level $\mathrm{\ell }=1,\dots ,d$, replace each consecutive pair of current nodes $(a_{\mathrm{\ell }-1,2j-1},a_{\mathrm{\ell }-1,2j})$, in increasing $j=1,\dots ,2^{d-\mathrm{\ell }}$, by $a_{\mathrm{\ell },j}=\mathrm{fl}(a_{\mathrm{\ell }-1,2j-1}+a_{\mathrm{\ell }-1,2j})\mathrm{.}$ If an operation fails, return range-failure with the first failing level $\mathrm{\ell }$ and node $j$ in this level-by-level, left-to-right order. For $d=0$, return $x_1$ with no additions. Otherwise return the single node $a_{d,1}$.

For termination and the distinction between the successful full-tree addition count and a failing initial segment, apply Finite arithmetic schedules terminate, with early failure counted separately to the prescribed levels. The operationwise model supplies each addition's finite value-or-failure return. Accuracy bounds require successful arithmetic.

For real floating-point a,b compute s=fl(a+b), z=fl(s-a), e=fl((a-(s-z))+(b-z)), with every operation separately rounded. Return the first range failure and operation index if an exact intermediate is nonzero subnormal, overflowed, or otherwise invalid; otherwise return (s,e). The fixed sequence terminates. Its general exact-transform claim is a separate theorem; this definition alone does not prove exact recovery.

For real floating-point a,b compute p=fl(ab), then e=fma(a,b,-p). Return range-failure if the product overflows, the exact FMA residual is nonzero subnormal or overflow, or another operationwise range condition fails; otherwise return (p,e).

The separate Finite arithmetic schedules terminate, with early failure counted separately applies to the ordered multiplication and FMA, each returning a finite result or immediate range failure. It gives at most two attempted rounded operations and exactly two on success. Exact recovery under representability of the residual as zero or a normal number is proved in TwoProduct with fused multiply-add returns the exact product residual.

For a finite list of $n\ge 0$ real floating-point inputs $x_1,\dots ,x_n$, initialize $s_0=c_0=0$. For $i=1,\dots ,n$ compute $y_i=\mathrm{fl}(x_i-c_{i-1})$, $s_i=\mathrm{fl}(s_{i-1}+y_i)$, $d_i=\mathrm{fl}(s_i-s_{i-1})$, and $c_i=\mathrm{fl}(d_i-y_i)$, in that order. Return $s_n$, or the first range-failure and its operation index.

The separate Finite arithmetic schedules terminate, with early failure counted separately applies to the four ordered rounded operations at each input index. It gives finite execution with at most $4n$ attempted rounded operations and exactly $4n$ on success, including zero for $n=0$. The successful error identity is proved in Kahan and Neumaier summation satisfy explicit finite data-dependent error bounds.

For a finite list of $n\ge 0$ real floating-point inputs $x_1,\dots ,x_n$, initialize $s_0=c_0=0$. For each $i=1,\dots ,n$, set $t_i=\mathrm{fl}(s_{i-1}+x_i)$. If $\mathrm{\mid }{s}_{i-1}\mathrm{\mid }\ge \mathrm{\mid }{x}_i\mathrm{\mid }$, set $q_i=\mathrm{fl}(s_{i-1}-t_i)$ and $a_i=\mathrm{fl}(q_i+x_i)$; otherwise set $q_i=\mathrm{fl}(x_i-t_i)$ and $a_i=\mathrm{fl}(q_i+s_{i-1})$. Then set $c_i=\mathrm{fl}(c_{i-1}+a_i)$ and $s_i=t_i$. Return $\mathrm{fl}(s_n+c_n)$, or the first range-failure and operation index.

The separate Finite arithmetic schedules terminate, with early failure counted separately applies to four rounded operations per input index, either branch, and the final addition. It gives finite execution with at most $4n+1$ attempted rounded operations and exactly $4n+1$ on success. The branch comparison is of initialized real values. Successful error properties are proved in Kahan and Neumaier summation satisfy explicit finite data-dependent error bounds.

Input a Hermitian positive-definite matrix $A\in {\mathbb{C}}^{n\times n}$, $n\ge 1$, whose real and imaginary components are in the floating-point system, using the convention of Hermitian positive-definite matrices and Cholesky factorisation A = LL* with positive diagonal; real matrices are included. Initialize every entry of $\widehat{L}$ to zero. Process columns $j=1,\dots ,n$ in increasing order. For the diagonal, initialize $s=0$ and, for $k=1,\dots ,j-1$ in increasing order, compute $r_k=\mathrm{fl}(\mathrm{\Re }({\widehat{l}}_{jk}{)}^2),t_k=\mathrm{fl}(\mathrm{\Im }({\widehat{l}}_{jk}{)}^2),q_k=\mathrm{fl}(r_k+t_k),s=\mathrm{fl}(s+q_k),$ in that order. Then compute $p_j=\mathrm{fl}(a_{jj}-s)$. If an operation fails, or if $p_j$ is not positive and normal, return failure with pivot $j$ and the first failing operation location. Otherwise compute ${\widehat{l}}_{jj}=\mathrm{fl}(\sqrt{p_j})$; if it is not positive and normal, return failure at pivot $j$.

For each $i=j+1,\dots ,n$ in increasing order, initialize the complex accumulator $c_{ij}=0$. For $k=1,\dots ,j-1$ in increasing order, compute the complex product $w_{ik}={\widehat{l}}_{ik}\overline{{\widehat{l}}_{jk}}$ using exactly the six real operations of the complex-product convention in Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add, then update $\mathrm{\Re }(c_{ij})\leftarrow \mathrm{fl}(\mathrm{\Re }(c_{ij})+\mathrm{\Re }(w_{ik})),\mathrm{\Im }(c_{ij})\leftarrow \mathrm{fl}(\mathrm{\Im }(c_{ij})+\mathrm{\Im }(w_{ik}))\mathrm{.}$ After this dot product, compute $r_{ij}=a_{ij}-c_{ij}$ by one rounded real subtraction for each component, then compute each component of ${\widehat{l}}_{ij}=r_{ij}\mathrm{/}{\widehat{l}}_{jj}$ by one rounded real division. The empty dot product is zero. On any failure return failure with pivot $j$ and its first failing operation location (stage, column, row, inner index, component, and operation where applicable), with the diagonal and row loops ordered as stated. If all columns finish, return success(\widehat L). By A matrix admits a Cholesky factorisation with positive diagonal exactly when it is Hermitian positive definite, and that factor is unique, exact positive definiteness ensures an exact factor exists, but does not ensure the rounded algorithm succeeds.

The finite column, row and inner-index loops unroll to an initialized arithmetic schedule. Each rounded operation returns a value or failure under the operationwise model, and each pivot test compares already computed values. The separate result Finite arithmetic schedules terminate, with early failure counted separately supplies termination with the specified completion or first-failure return; no arithmetic-success or backward-error conclusion follows from termination alone.

Let $\mathbb{F}\in \{\mathbb{R},\mathbb{C}\}$. Input a real active column $x\in {\mathbb{R}}^r\subseteq {\mathbb{F}}^r$ and trailing block $B\in {\mathbb{F}}^{r\times q}$ with entries in the floating-point system, $r\ge 1$, $q\ge 0$. If every component of $x$ is zero, return the identity action and a copy of $B$, with the active output column set to exact zeros, so that $C_{+}=[0e_{1}B]$; return $(v=0,\widehat{\tau }=0,\widehat{\alpha }=0,C_{+})$ marked identity. Otherwise compute $s={\max }_{1\le i\le r}\mathrm{\mid }{x}_i\mathrm{\mid }>0$, and then, in increasing $i$, compute $z_i=\mathrm{fl}(x_i\mathrm{/}s)$. Initialize $t=0$ and, in increasing $i$, compute $a_i=\mathrm{fl}(z_i{z}_i)$ followed by $t=\mathrm{fl}(t+a_i)$. Compute $h=\mathrm{fl}(\sqrt{t})$ and then $\widehat{r}=\mathrm{fl}(sh)$. Set $\sigma =1$ if $x_1=0$ and $\sigma =\mathrm{sign}(x_1)$ otherwise, and set $\widehat{\alpha }=-\sigma \widehat{r}$ using the exact sign change. Compute $v_1=\mathrm{fl}(x_1-\widehat{\alpha })$, then copy $v_i=x_i$ exactly for $i=2,\dots ,r$.

Initialize $w=0$ and, in increasing $i$, compute $a_i=\mathrm{fl}(v_i{v}_i)$ followed by $w=\mathrm{fl}(w+a_i)$. Then compute $\widehat{\tau }=\mathrm{fl}(2\mathrm{/}w)$. For each column $j=1,\dots ,q$ in increasing order, initialize $d_j=0$ and, in increasing $i$, compute $a_{ij}=\mathrm{fl}(v_i{B}_{ij})$ followed by $d_j=\mathrm{fl}(d_j+a_{ij})$. After all $d_j$ have been computed, in increasing $i$ compute $p_i=\mathrm{fl}(\widehat{\tau }{v}_i)$. Then, in increasing $j$ and within each column in increasing $i$, compute $U_{ij}=\mathrm{fl}(p_i{d}_j)$ and ${\widehat{B}}_{ij}=\mathrm{fl}(B_{ij}-U_{ij})$ in that order. The complete active output is $C_{+}=[\widehat{\alpha }{e}_1\widehat{B}]$, with exact zeros written below its first diagonal entry. If any elementary operation fails, return the first range failure with its stage and indices in the order above. Also return failure if $t$, $\widehat{r}$, or $w$ is not positive normal, or any intermediate is nonfinite. Otherwise return $(v,\widehat{\tau },\widehat{\alpha },C_{+})$ marked reflector. For a separate block $Y\in {\mathbb{F}}^{r\times q^{\mathrm{\prime }}}$, $q^{\mathrm{\prime }}\ge 0$, apply the stored reflector by computing its sequential column dots in increasing column order, then the vector $p_i=\mathrm{fl}(\widehat{\tau }{v}_i)$ in increasing $i$, and then each rounded product and subtraction by increasing column and row; omit the target-column overwrite. Store reflector factors with their identity/reflector marker to represent the exact orthogonal product; a rounded explicit dense product is a separate output with its own arithmetic error. For a nonzero stored $v$, the associated exact reflector is $H(v)=I-2v{v}^T\mathrm{/}(v^{T}v)$, using Householder reflectors in real or complex inner-product spaces and the QR convention of Full, reduced, and column-pivoted computational QR factorisations. Since $v$ is real, this matrix acts on ${\mathbb{F}}^r$ as an orthogonal reflector over $\mathbb{R}$ and a unitary reflector over $\mathbb{C}$. This definition does not assert a local backward-error bound.

On every admitted finite input of Scaled Householder reflector formation and application, the formation and application procedure terminates with its specified identity, reflector or first-failure return. Each specified separate-block application also terminates with its output or first failure. The zero-column branch returns the identity action and the copied block without evaluating a reflector denominator. In a successful nonzero-column branch the stored vector is nonzero, so its associated exact real reflector is orthogonal and also acts unitarily on complex blocks. This statement asserts neither unconditional arithmetic success nor a backward-error bound for a successful call.

The reference entries are inherited from the analyzed procedure as method context. This result extracts and composes its explicit finite schedule and representative algebra; no verbatim external theorem or additional external reading is claimed.

This is a local extraction of finite-loop and representative arguments from the existing procedure carrier. Its error estimates remain in lem-scaled-householder-step-meets-local-backward-error-premise under that lemma's distinct hypotheses; no external full-text reading is claimed.

Input nonsingular real $A\in {\mathbb{R}}^{n\times n}$, $n\ge 1$, exact data $b\in {\mathbb{R}}^n$, and an initial iterate $x_0$ stored at solve/update precision $u_w$. Also input stored triangular factors $\widehat{L},\widehat{U}$ formed at factorization precision $u_f$, where $\widehat{L}$ is lower triangular, $\widehat{U}$ is upper triangular, both have nonzero diagonals, and $\widehat{A}:=\widehat{L}\widehat{U}$ is the exact real product; thus $\widehat{A}$ is nonsingular. Input residual precision $u_r$, solve/update precision $u_w$, finite update budget $N\ge 0$, tolerance $\tau \ge 0$, a stagnation window $K_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}\ge 1$, nonnegative constants ${\eta }_r,{\xi }_r,{\eta }_w,{\xi }_w$, and a proposed contraction upper bound $0\le q_{\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{t}}<1$. Let $x_{\ast }=A^{-1}b$. Use the operationwise model of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add at each stated precision. Entries of the exact data $A,b$ are available as real operands to the residual operations; the stored factors, iterates, residuals, and corrections are used as their exact represented values in those operations.

Certificates have the following fixed format: one real matrix $B\in {\mathbb{R}}^{n\times n}$ before iteration, nonnegative real scalars ${\rho }_k$ for $0\le k\le N$, and nonnegative real scalars ${\omega }_k$ for $0\le k<N$. Each entry may instead be marked missing. All certificate checks use exact real addition, subtraction, multiplication, division by a positive number, nonnegative square roots, and comparisons; they are not rounded operations at $u_r$ or $u_w$. Evaluate finite sums and products in increasing index order. In particular, evaluate $\mathrm{\parallel }v{\mathrm{\parallel }}_2=\sqrt{{\sum }_i{v}_i^2}$ and $\mathrm{\parallel }M{\mathrm{\parallel }}_F=\sqrt{{\sum }_{i,j}{M}_{ij}^2}$ by these exact operations. This is an ideal exact-real certificate verifier, not a claim that arbitrary real inputs admit a finite-bit implementation.

Before iterating, check triangularity and nonzero diagonals of the factors; on failure return inconclusive with $x_0$. A missing or incorrectly shaped $B$ returns inconclusive with $x_0$. Form $\widehat{A}=\widehat{L}\widehat{U}$ exactly and check $B\widehat{A}=I$ entrywise. If this fails, return inconclusive with $x_0$. Set $a_{+}=\mathrm{\parallel }A{\mathrm{\parallel }}_F,a_{-}=a_{+}\mathrm{/}\sqrt{n},h=\mathrm{\parallel }B{\mathrm{\parallel }}_F,d=\mathrm{\parallel }I-BA{\mathrm{\parallel }}_F\mathrm{.}$ Accept the initial certificate precisely when $d+h((1+{\eta }_w){\eta }_r+{\eta }_w)a_{+}\le q_{\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{t}}<1;$ otherwise return inconclusive with $x_0$. Here $B={\widehat{A}}^{-1}$, and the inequalities $a_{-}\le \mathrm{\parallel }A{\mathrm{\parallel }}_2\le a_{+}$ and $\mathrm{\parallel }M{\mathrm{\parallel }}_2\le \mathrm{\parallel }M{\mathrm{\parallel }}_F$ from Spectral and Frobenius norms are unitarily invariant, are given by singular values, and satisfy the sharp rank comparison make this a sufficient certificate for $q=\mathrm{\parallel }I-{\widehat{A}}^{-1}A{\mathrm{\parallel }}_2+\mathrm{\parallel }{\widehat{A}}^{-1}{\mathrm{\parallel }}_2((1+{\eta }_w){\eta }_r+{\eta }_w)\mathrm{\parallel }A{\mathrm{\parallel }}_2\le q_{\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{t}}\mathrm{.}$ These conservative tests may reject data for which sharper norm certificates exist.

Before entering the iteration loop, initialize the consecutive-equality counter $c_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}=0$ and mark the preceding residual as undefined. Retain this counter and preceding residual across iterations; do not reinitialize them inside the loop. The counter records the number of consecutive equalities between adjacent successfully certified computed residual vectors.

For each $k=0,\dots ,N$, compute the residual by rows in increasing $i$: initialize $s_i=0$, then for $j=1,\dots ,n$ in increasing order compute $p_{ij}={\mathrm{fl}}_{u_r}(A_{ij}{x}_{k,j})$ and $s_i={\mathrm{fl}}_{u_r}(s_i+p_{ij})$; finally set ${\widehat{r}}_{k,i}={\mathrm{fl}}_{u_r}(b_i-s_i)$. Any operationwise range failure returns inconclusive with the current iterate and the first failing row, column, and operation. Define the exact residual defect ${\delta }_{r,k}={\widehat{r}}_k-(b-A{x}_k)$, so ${\widehat{r}}_k=A(x_{\ast }-x_k)+{\delta }_{r,k}$. Check that ${\rho }_k$ is present and nonnegative and that ${\sum }_i{\delta }_{r,k,i}^2\le {\rho }_k^2,{\rho }_k\le {\eta }_r(a_{-}\mathrm{\parallel }{x}_k{\mathrm{\parallel }}_2+\mathrm{\parallel }b{\mathrm{\parallel }}_2)+{\xi }_r\mathrm{.}$ A missing scalar or failed comparison returns inconclusive with the current iterate. These exact tests imply $\mathrm{\parallel }{\delta }_{r,k}{\mathrm{\parallel }}_2\le {\eta }_r(\mathrm{\parallel }A{\mathrm{\parallel }}_2\mathrm{\parallel }{x}_k{\mathrm{\parallel }}_2+\mathrm{\parallel }b{\mathrm{\parallel }}_2)+{\xi }_r$. Use the retained consecutive-equality counter and preceding residual. Use the exact Euclidean norm of the stored residual vector for the status comparison; this is a comparison of the already computed vector, not an additional rounded arithmetic operation. If $\mathrm{\parallel }{\widehat{r}}_k{\mathrm{\parallel }}_2\le \tau$, return residual-tolerance with $x_k$. If a preceding residual exists, compare its components exactly with those of ${\widehat{r}}_k$: if all are equal, increment $c_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}$ by one; otherwise reset it to zero. If $c_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}=K_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}$, return successful-stagnation with $x_k$. This status certifies only that the last $K_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{g}}$ successful updates produced identical computed residual vectors, with all range and contraction certificates accepted; it certifies neither a small exact error nor convergence. Before any further update, store ${\widehat{r}}_k$ as the preceding residual.

If $k=N$ and neither status above was returned, return budget-exhausted with $x_k$. Otherwise solve $\widehat{L}{y}_k={\widehat{r}}_k$ by forward substitution, then $\widehat{U}{\widehat{z}}_k=y_k$ by backward substitution, using precision $u_w$. For forward substitution, visit rows $i=1,\dots ,n$ in increasing order, initialize $s_i=0$, and for $j=1,\dots ,i-1$ in increasing order compute $p_{ij}={\mathrm{fl}}_{u_w}({\widehat{L}}_{ij}{y}_{k,j})$ and $s_i={\mathrm{fl}}_{u_w}(s_i+p_{ij})$, then $t_i={\mathrm{fl}}_{u_w}({\widehat{r}}_{k,i}-s_i)$ and $y_{k,i}={\mathrm{fl}}_{u_w}(t_i\mathrm{/}{\widehat{L}}_{ii})$. For backward substitution, visit rows $i=n,\dots ,1$ in decreasing order, initialize $s_i=0$, and for $j=i+1,\dots ,n$ in increasing order compute $p_{ij}={\mathrm{fl}}_{u_w}({\widehat{U}}_{ij}{\widehat{z}}_{k,j})$ and $s_i={\mathrm{fl}}_{u_w}(s_i+p_{ij})$, then $t_i={\mathrm{fl}}_{u_w}(y_{k,i}-s_i)$ and ${\widehat{z}}_{k,i}={\mathrm{fl}}_{u_w}(t_i\mathrm{/}{\widehat{U}}_{ii})$. Update components in increasing order by $x_{k+1,i}={\mathrm{fl}}_{u_w}(x_{k,i}+{\widehat{z}}_{k,i})$. A failed operation returns inconclusive with the current iterate and first failing operation location.

After a successful update, define the exact residual-equivalent update defect ${\delta }_{w,k}:=\widehat{A}(x_{k+1}-x_k)-{\widehat{r}}_k,$ so $x_{k+1}=x_k+{\widehat{A}}^{-1}({\widehat{r}}_k+{\delta }_{w,k})$. Check that ${\omega }_k$ is present and nonnegative and that ${\sum }_i{\delta }_{w,k,i}^2\le {\omega }_k^2$ and ${\omega }_k\le {\eta }_w\mathrm{\parallel }{\widehat{r}}_k{\mathrm{\parallel }}_2+{\xi }_w$; otherwise return inconclusive with $x_{k+1}$. These tests imply $\mathrm{\parallel }{\delta }_{w,k}{\mathrm{\parallel }}_2\le {\eta }_w\mathrm{\parallel }{\widehat{r}}_k{\mathrm{\parallel }}_2+{\xi }_w$. Continue until one of the stated statuses is returned. All loops are finite. The procedure does not claim exact termination or unconditional convergence.

Let $x_1,\dots ,x_n\in \mathbb{R}$, $n\ge 1$, be successfully processed by Balanced pairwise summation with range-failure status in the operationwise model. Put $d=\lceil {\log }_{2}n\rceil$, assume $du<1$, and set ${\gamma }_d=du\mathrm{/}(1-du)$. Then $\mid \widehat{s}-{\sum }_{i=1}^n{x}_i\mid \le {\gamma }_d{\sum }_{i=1}^n\mathrm{\mid }{x}_i\mathrm{\mid }\mathrm{.}$ For $n=1$, $d=0$ and the error is zero.

For an integer $p\ge 2$, let ${\mathbb{B}}_p=\{0\}\cup \{\pm m{2}^k:m,k\in \mathbb{Z},2^{p-1}\le m<2^p\}$, with unbounded exponent. Let $R$ round a real number to its nearest member of ${\mathbb{B}}_p$, breaking ties by the even normalized significand; set $R(0)=0$. For $x\mathrm{\ne }0$ put $E(x)$ be the unique integer with $2^{E(x)}\le \mathrm{\mid }x\mathrm{\mid }<2^{E(x)+1}$, and set $q(x)=2^{E(x)-p+1}$. The following hold.

• $R$ is monotone, odd, and fixes ${\mathbb{B}}_p$. For $R(t)\mathrm{\ne }0$, $\mathrm{\mid }R(t)-t\mathrm{\mid }\le q(R(t))\mathrm{/}2$.
• If $x,y\in {\mathbb{B}}_p$, then $\rho =R(x+y)-(x+y)$ belongs to ${\mathbb{B}}_p$ and $\mathrm{\mid }\rho \mathrm{\mid }\le \min (\mathrm{\mid }x\mathrm{\mid },\mathrm{\mid }y\mathrm{\mid })$. If $\mathrm{\mid }R(x+y)\mathrm{\mid }<\min (\mathrm{\mid }x\mathrm{\mid },\mathrm{\mid }y\mathrm{\mid })$, the sum is exact.
• If $x,y\in {\mathbb{B}}_p$ have the same sign and $\mathrm{\mid }x\mathrm{\mid }\mathrm{/}2\le \mathrm{\mid }y\mathrm{\mid }\le 2\mathrm{\mid }x\mathrm{\mid }$, then $x-y\in {\mathbb{B}}_p$ (Sterbenz exact subtraction).
• If $A,B\in {\mathbb{B}}_p$ and $\mathrm{\mid }A\mathrm{\mid }\ge \mathrm{\mid }B\mathrm{\mid }$, then $S=R(A+B)$, $Z=R(S-A)$, $T=R(B-Z)$ satisfy $Z=S-A$ and $T=A+B-S$ exactly.

Use a binary floating-point format of precision $p\ge 2$: its normal positive values are $m{2}^{j-p+1}$ for integers $2^{p-1}\le m<2^p$ and $e_{\min }\le j\le e_{\max }$, with their negatives and zero. Any additional representable values below the normal range must belong to the unbounded binary set ${\mathbb{B}}_p$ of Binary rounding grids, exact subtraction, and ordered TwoSum (as ordinary binary subnormals do). Round each operation to nearest, ties to even, with the no-underflow success rule of Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add.

For any representable real inputs $a,b$ for which TwoSum error-free transformation succeeds, its output satisfies $s=\mathrm{fl}(a+b),e=(a+b)-s,s+e=a+b$ exactly in real arithmetic. No ordering or sign condition on the inputs is required; zeros, exact cancellation, and all successful calls at either normal-range boundary are included.

The binary digit hypothesis supplies exactness information absent from the relative-error model alone. Shewchuk's Theorem 7 states the binary result for $p\ge 3$; step 3.1 uses the sharper bound on $s-a$ directly to include $p=2$ as well. The theorem retains unrestricted inputs and does not assume that the residual operations are exact: their exactness is proved above.

Let $a,b$ be real floating-point numbers successfully processed by TwoProduct product and residual computation using fused multiply-add, and assume the exact residual $ab-p$ is representable as zero or a normal number. Then the returned pair satisfies $p+e=ab$ exactly.

Let the input length be $n\in {\mathbb{N}}_0$. If $n=0$, both algorithms return exactly zero, the sum of the empty input, and their absolute errors are zero. Kahan executes no loop operations and returns $s_0=0$; Neumaier executes no loop operations and its final addition returns $\mathrm{fl}(s_0+c_0)=0$, with final rounding error ${ϵ}^f=0$. The Kahan identity and bound below apply when $n\ge 1$; no $y_0$, $s_{-1}$, or ${ϵ}_0^s$ is required in the empty branch.

For $n\ge 1$ successfully processed real inputs $x_1,\dots ,x_n$, let ${ϵ}_i^y,{ϵ}_i^s,{ϵ}_i^d,{ϵ}_i^c$ be the signed errors in the four Kahan operations. Then $s_n-{\sum }_{i=1}^n{x}_i={\sum }_{i=1}^n{ϵ}_i^y+{ϵ}_n^s-{\sum }_{i=1}^{n-1}({ϵ}_i^d+{ϵ}_i^c),$ hence $\mathrm{\mid }{s}_n-{\sum }_i{x}_i\mathrm{\mid }\le u{\sum }_i\mathrm{\mid }{x}_i-c_{i-1}\mathrm{\mid }+u\mathrm{\mid }{s}_{n-1}+y_n\mathrm{\mid }+u{\sum }_{i<n}(\mathrm{\mid }{s}_i-s_{i-1}\mathrm{\mid }+\mathrm{\mid }{d}_i-y_i\mathrm{\mid })\mathrm{.}$ For Neumaier, let ${ϵ}_i^t,{ϵ}_i^q,{ϵ}_i^a,{ϵ}_i^c$ be the signed errors at the actually executed operations and ${ϵ}^f$ the final rounding error. Then $\widehat{s}-{\sum }_i{x}_i={\sum }_i({ϵ}_i^q+{ϵ}_i^a+{ϵ}_i^c)+{ϵ}^f\mathrm{.}$ The leading-sum error ${ϵ}_i^t$ cancels from this exact identity; its nonnegative allowance may still be retained in the conservative bound below. For each $i$, define $P_i$ by the branch-specific expression below. Then $\mathrm{\mid }\widehat{s}-{\sum }_i{x}_i\mathrm{\mid }\le u{\sum }_i{P}_i+u\mathrm{\mid }{s}_n+c_n\mathrm{\mid },$ where the final term accounts for the last rounded addition. These are finite data-dependent bounds, not uniform relative-error claims.

$P_i=\mathrm{\mid }{s}_{i-1}\mathrm{\mid }+\mathrm{\mid }{x}_i\mathrm{\mid }+\{\begin{array}{ll}\mathrm{\mid }{s}_{i-1}\mathrm{\mid }+\mathrm{\mid }{t}_i\mathrm{\mid },& \mathrm{\mid }{s}_{i-1}\mathrm{\mid }\ge \mathrm{\mid }{x}_i\mathrm{\mid },\\ \mathrm{\mid }{x}_i\mathrm{\mid }+\mathrm{\mid }{t}_i\mathrm{\mid },& \mathrm{\mid }{s}_{i-1}\mathrm{\mid }<\mathrm{\mid }{x}_i\mathrm{\mid }\end{array}+\{\begin{array}{ll}\mathrm{\mid }{q}_i\mathrm{\mid }+\mathrm{\mid }{x}_i\mathrm{\mid },& \mathrm{\mid }{s}_{i-1}\mathrm{\mid }\ge \mathrm{\mid }{x}_i\mathrm{\mid },\\ \mathrm{\mid }{q}_i\mathrm{\mid }+\mathrm{\mid }{s}_{i-1}\mathrm{\mid },& \mathrm{\mid }{s}_{i-1}\mathrm{\mid }<\mathrm{\mid }{x}_i\mathrm{\mid }\end{array}+\mathrm{\mid }{c}_{i-1}\mathrm{\mid }+\mathrm{\mid }{a}_i\mathrm{\mid }\mathrm{.}$

Let $n\ge 1$ and let $A=A^{\ast }$ be positive definite in ${\mathbb{C}}^{n\times n}$. Suppose the rounded Cholesky algorithm succeeds with $\widehat{L}$. Put $m=20n+20$ and $g={\gamma }_m=\frac{mu}{1-mu}$; assume $mu<1\mathrm{/}5$. For $E=\widehat{L}{\widehat{L}}^{\ast }-A$, which is Hermitian, and with $\mathrm{\mid }M\mathrm{\mid }$ denoting entrywise absolute value for matrices, $\mathrm{\mid }E\mathrm{\mid }\le 100g\mathrm{\mid }\widehat{L}\mathrm{\mid }\mathrm{\mid }\widehat{L}{\mathrm{\mid }}^{\ast }\text{entrywise},\mathrm{\parallel }E{\mathrm{\parallel }}_F\le 100g\mathrm{\parallel }\widehat{L}{\mathrm{\parallel }}_F^2,\mathrm{\parallel }E{\mathrm{\parallel }}_2\le 100g\mathrm{\parallel }\widehat{L}{\mathrm{\parallel }}_F^2\mathrm{.}$ Thus $c_n(u)u=\frac{100mu}{1-mu}$. If also $100\sqrt{n}g\le 1\mathrm{/}2$, then $\mathrm{\parallel }E{\mathrm{\parallel }}_F\le 200\sqrt{n}g\mathrm{\parallel }A{\mathrm{\parallel }}_F,\mathrm{\parallel }E{\mathrm{\parallel }}_2\le 200\sqrt{n}g\mathrm{\parallel }A{\mathrm{\parallel }}_F\mathrm{.}$ These bounds apply on success only; positive definiteness does not guarantee rounded success.