Numerical Stability

Why Numerical Stability Matters

A sensor logs one enormous baseline, two small adjustments, and the negative of the baseline. The exact total of the stored numbers $[2^{53},1,1,-2^{53}]$ is $2$. A straightforward floating-point loop returns $0$: each adjustment disappears before the baseline is cancelled. The program can follow the formula perfectly and still give an answer that is useless for the application.

Numerical stability studies how arithmetic errors propagate through an algorithm. This lesson focuses on summation because every dot product, matrix operation, statistical accumulator, and many simulations contain a reduction to a total. You will trace rounding one operation at a time, implement Kahan and Neumaier compensation, compare a balanced reduction, and use partial pivoting and an independently evaluated residual to check a small linear solve.

The same habits apply beyond sums: identify the sensitive operation, choose a better formulation, specify the arithmetic assumptions, and test against a reference with more trustworthy arithmetic. An algorithm cannot recover information already lost when the inputs were measured or converted to floating point.

Core Concepts and Terminology

The input to a summation routine is an ordered finite sequence of stored floating-point values $x_1,\ldots,x_n$. Its mathematical target is the exact real sum $S=\sum_{i=1}^{n}x_i$ of those stored values; its output $\widehat S$ is another floating-point value. This target must be distinguished from the sum of the intended decimal quantities before conversion. An empty sequence has sum $0$.

Assume IEEE 754 binary64, round to nearest with ties to even, and separately rounded operations. A normal binary64 number has $53$ significant binary digits, including the implicit leading digit. The gap above $1$ is $\epsilon=2^{-52}$, while the unit roundoff used here is $u=2^{-53}$; authors sometimes use “machine epsilon” for either quantity, so check the convention. Spacing grows with magnitude. The next representable number above $M=2^{53}=9007199254740992$ is $M+2$, and the exact midpoint $M+1$ rounds to $M$.

For a basic operation whose exact result stays in the normal range, with no overflow, the usual model is

$$\operatorname{fl}(a\mathbin{\circ}b)=(a\mathbin{\circ}b)(1+\delta),\qquad |\delta|\le u.$$

Here $\circ$ denotes addition, subtraction, multiplication, or division, with a nonzero divisor. This is a model for rounding one operation on the stored operands; it does not say the operands are accurate measurements. Subnormal results need an absolute-error treatment near the minimum spacing $2^{-1074}$, rather than this uniform relative bound. NaN, infinity, overflow, and an environment that flushes subnormals to zero require separate handling.

Two distinctions prevent common misunderstandings. First, floating-point addition is not associative: $\operatorname{fl}(\operatorname{fl}(M+1)-M)=0$, whereas $\operatorname{fl}(M+\operatorname{fl}(1-M))=1$. In the second expression, $1-M$ is exactly representable because spacing just below $M$ is $1$.

Second, cancellation need not create a rounding error in the subtraction itself. Two nearby stored numbers can subtract exactly, yet their difference may be a poor approximation to the difference of the original unrounded inputs. For example, first storing $M+1$ as $M$, then subtracting $M$, produces an exact subtraction yielding $0$ instead of the intended $1$. Compensating that subtraction cannot reconstruct the discarded input bit.

For summation, define the componentwise relative condition number, when $S\ne0$, by

$$\kappa_{\mathrm{sum}}=\frac{\sum_i |x_i|}{|S|}.$$

If each input changes by at most $\eta|x_i|$, the exact total changes by at most $\eta\sum_i|x_i|$, giving relative change at most $\eta\kappa_{\mathrm{sum}}$. Terms with one sign have condition number $1$; strong cancellation can make it enormous. The worked sequence has $\kappa_{\mathrm{sum}}=M+1$. If $S=0$, relative output error is undefined, so use an absolute measure instead; if all terms vanish, the sum itself is trivial.

How Numerical Stability Works, Step by Step

Start with the naive algorithm: set $s=0$ and repeatedly replace it by $\operatorname{fl}(s+x_i)$. Every rounded total becomes the input to the next addition. Once a small term vanishes, this algorithm keeps no record of it.

Kahan summation maintains $s$ and $c$, both initially zero. For each stored input $x$, perform these four assignments in this exact order, with every displayed arithmetic operation rounded separately:

1. $y=\operatorname{fl}(x-c)$: adjust the new term using the previous correction.
2. $t=\operatorname{fl}(s+y)$: form the candidate total.
3. $c=\operatorname{fl}(\operatorname{fl}(t-s)-y)$: estimate the signed error in this addition.
4. $s=t$: commit the total only after the old $s$ has been used to compute $c$.

Return $s$ after the final input. In this version, the correction is subtracted from the next input; another formulation may define its correction with the opposite sign. Do not mix the formulas. Returning $s+c$ here is also a different algorithm, not an extra improvement to this convention.

Why do apparently redundant operations help? When $s=M$ and $y=1$, the rounded candidate is $t=M$. The calculation $(t-s)-y$ becomes $0-1=-1$, recording that the new total fell short of the intended addition by $1$. The next input $1$ is adjusted to $1-(-1)=2$, large enough to survive addition to $M$.

flowchart TD
    A["Set s = 0 and c = 0"] --> B{"Another input?"}
    B -->|Yes| C["Read finite x; compute y = x - c"]
    C --> D["Compute t = s + y"]
    D --> E["Compute c = (t - s) - y using old s"]
    E --> F["Set s = t"]
    F --> B
    B -->|No| G["Return s"]

The implementations also reject nonfinite inputs and nonfinite intermediate results. Their arithmetic remains finite for the worked examples. This check catches overflow instead of silently continuing with infinity, but it does not detect every underflow or every loss of accuracy.

Neumaier summation uses a related state with a different sign convention. Set $t=\operatorname{fl}(s+x)$ and extract the addition's omitted part using the larger-magnitude operand first:

$$e=\begin{cases}\operatorname{fl}(\operatorname{fl}(s-t)+x),&|s|\ge|x|,\\\operatorname{fl}(\operatorname{fl}(x-t)+s),&|s|\lt|x|.\end{cases}$$

Then set $c=\operatorname{fl}(c+e)$ and $s=t$; return $\operatorname{fl}(s+c)$ at the end. Under the usual error-free addition-transform assumptions, including no overflow and ordinary gradual underflow, the chosen branch recovers an $e$ satisfying the exact real identity $s_{\mathrm{old}}+x=t+e$. Accumulating all such corrections and combining them with $s$ still incurs rounding, so the complete sum is not universally exact.

Balanced pairwise summation offers a third strategy: split the index interval in half, sum both halves recursively, then add their totals. Its advantage comes from reducing the number of additions on an input-to-output path, not from retaining a correction. Fixing this tree matters when reproducible results are required.

Worked Example with a Full Trace

Use $M=2^{53}$ and the ordered values $[M,1,1,-M]$. All four inputs and the exact answer $2$ are representable. Thus this example isolates algorithmic rounding error rather than decimal-conversion error.

The naive totals after each input are $M,M,M,0$. Both additions of $1$ encounter the same halfway case and round down to $M$. The final subtraction is exact; it merely reveals that the two earlier contributions were lost.

For Kahan, the following table records every arithmetic intermediate. Each row uses the previous row's $s$ and $c$; the final column is the newly committed $s$.

On the last row, $2-(M+2)=-M$ is exact. The algorithm returns $2$, so its absolute error on these stored inputs is zero. This successful trace is an example, not a proof that all Kahan sums are exact.

For comparison, Neumaier's successive states $(s,c)$ are $(M,0)$, $(M,1)$, $(M,2)$, and $(0,2)$, followed by the final addition $s+c=2$. It stores positive omitted pieces in $c$, whereas this Kahan version carries a negative correction into the next input.

The specified pairwise tree produces $\operatorname{fl}(M+1)=M$ in its left half and $\operatorname{fl}(1-M)=1-M$ in its right half. Adding those two results gives $1$. It improves the naive answer here but still misses the exact total. A better worst-case error bound does not imply a better answer on every individual input.

A second application concerns the system

$$A=\begin{bmatrix}\varepsilon&1\\1&1\end{bmatrix},\qquad b=\begin{bmatrix}1\\2\end{bmatrix},\qquad\varepsilon=10^{-20}.$$

Here $\varepsilon$ denotes the stored approximation to $10^{-20}$; it is unrelated to the spacing $\epsilon=2^{-52}$ defined earlier. Without pivoting, the first multiplier is approximately $10^{20}$. Both $1-10^{20}$ and $2-10^{20}$ round to $-10^{20}$. Back substitution then gives $\widehat x=[0,1]$, whose residual infinity norm is $1$.

Partial pivoting first swaps the rows, using a pivot of $1$ and a multiplier of about $10^{-20}$. The second row becomes $[0,1\mid1]$ after rounding, and back substitution gives $\widehat x=[1,1]$. Its exact residual on the stored system is $[-\varepsilon,0]$, so the infinity norm is about $10^{-20}$. The exact solution is $[1/(1-\varepsilon),(1-2\varepsilon)/(1-\varepsilon)]$, whose components round to $1$ in binary64. The matrix is well conditioned at this scale: the bad unpivoted answer comes from the algorithm, not from an inherently sensitive system.

Implementation in Python, C++, and Java

Each panel is a complete program. The four summation routines use identical operation ordering, and the naive routine is an explicit loop so that a language's built-in summation strategy cannot change the comparison. The balanced tree splits half-open index ranges at their midpoint and avoids copying slices. No summation routine sorts or mutates its input.

The linear solver uses plain nested arrays, copies the augmented system, selects the largest absolute entry in each active column, and swaps whole rows. It rejects a pivot if its magnitude divided by the largest absolute entry in the original matrix is at most pivot_tol. That guard is an operational scale test, not a proof of singularity or a condition-number estimate; it can reject a nonsingular system with badly mismatched scales. Row/column equilibration or a library solver is preferable for such systems.

Back substitution uses Neumaier accumulation for its already-rounded products. residual_inf independently evaluates $b-A\widehat x$ from the original inputs, including $b_i$ and the negative products in one compensated row sum. It improves accumulation, but product rounding remains. For this particular sample it retains the tiny $-\varepsilon$ term that a straightforward dot product followed by subtraction can hide. The printed residual check uses absolute tolerance $10^{-12}$ for this fixed, order-one example; it is not a universal tolerance.

Use Python 3 with binary64 floats, C++ with binary64 double, or Java 17 or later. Compile C++ with g++ -std=c++14 -O2 -fno-fast-math -ffp-contract=off main.cpp -o stability; run Java with javac Main.java followed by java Main. Preserve parentheses, separate rounding, and the rounding mode. Reassociation or unintended contraction can invalidate a compensation analysis and can change the demonstration.

import math
import sys


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


def naive_sum(values):
    s = 0.0
    for value in values:
        s = finite(s + finite(value))
    return s


def kahan_sum(values):
    s, c = 0.0, 0.0
    for value in values:
        y = finite(finite(value) - c)
        t = finite(s + y)
        c = finite((t - s) - y)
        s = t
    return s


def neumaier_sum(values):
    s, c = 0.0, 0.0
    for value in values:
        value = finite(value)
        t = finite(s + value)
        if abs(s) >= abs(value):
            error = finite((s - t) + value)
        else:
            error = finite((value - t) + s)
        c = finite(c + error)
        s = t
    return finite(s + c)


def pairwise_sum(values):
    def visit(lo, hi):
        if lo == hi:
            return 0.0
        if hi - lo == 1:
            return finite(values[lo])
        mid = lo + (hi - lo) // 2
        return finite(visit(lo, mid) + visit(mid, hi))
    return visit(0, len(values))


def solve_pivoted(a, b, pivot_tol=1e-14):
    n = len(b)
    if n == 0 or len(a) != n or any(len(row) != n for row in a):
        raise ValueError("expected a nonempty square system")
    if not math.isfinite(pivot_tol) or not 0.0 < pivot_tol < 1.0:
        raise ValueError("invalid pivot tolerance")
    rows = [[finite(v) for v in row] + [finite(rhs)]
            for row, rhs in zip(a, b)]
    scale = max(abs(rows[i][j]) for i in range(n) for j in range(n))
    if scale == 0.0:
        raise ValueError("zero matrix")
    for j in range(n):
        p = max(range(j, n), key=lambda i: abs(rows[i][j]))
        if abs(rows[p][j]) / scale <= pivot_tol:
            raise ValueError("pivot too small for this scale")
        rows[j], rows[p] = rows[p], rows[j]
        for i in range(j + 1, n):
            factor = finite(rows[i][j] / rows[j][j])
            rows[i][j] = 0.0
            for k in range(j + 1, n + 1):
                rows[i][k] = finite(rows[i][k] - factor * rows[j][k])
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        terms = [finite(rows[i][j] * x[j]) for j in range(i + 1, n)]
        rhs = finite(rows[i][n] - neumaier_sum(terms))
        x[i] = finite(rhs / rows[i][i])
    return x


def residual_inf(a, b, x):
    n = len(b)
    if n == 0 or len(a) != n or len(x) != n or any(len(row) != n for row in a):
        raise ValueError("incompatible residual dimensions")
    result = 0.0
    for row, rhs in zip(a, b):
        terms = [finite(rhs)]
        terms.extend(finite(-finite(v) * finite(z)) for v, z in zip(row, x))
        result = max(result, abs(neumaier_sum(terms)))
    return result


def main():
    if sys.float_info.radix != 2 or sys.float_info.mant_dig != 53:
        raise RuntimeError("this example requires binary64")
    m = 9007199254740992.0
    values = [m, 1.0, 1.0, -m]
    cancellation = [1e16, 1.0, -1e16]
    results = [naive_sum(values), kahan_sum(values), neumaier_sum(values),
               pairwise_sum(values), kahan_sum(cancellation), neumaier_sum(cancellation)]
    labels = ["trace naive", "trace kahan", "trace neumaier", "trace pairwise",
              "cancellation kahan", "cancellation neumaier"]
    for label, value in zip(labels, results):
        print(f"{label} = {value:.0f}")
    print("summation checks = " + ("PASS" if results == [0., 2., 2., 1., 0., 1.] else "FAIL"))
    a = [[1e-20, 1.0], [1.0, 1.0]]
    b = [1.0, 2.0]
    x = solve_pivoted(a, b)
    print(f"pivoted x = [{x[0]:.6f}, {x[1]:.6f}]")
    print("residual check = " + ("PASS" if residual_inf(a, b, x) <= 1e-12 else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
using Vector = std::vector<double>;
using Matrix = std::vector<Vector>;

double finite_value(double value) {
    if (!std::isfinite(value)) throw std::invalid_argument("nonfinite input or intermediate result");
    return value;
}

double naive_sum(const Vector& values) {
    double s = 0.0;
    for (double v : values) s = finite_value(s + finite_value(v));
    return s;
}

double kahan_sum(const Vector& values) {
    double s = 0.0, c = 0.0;
    for (double v : values) {
        double y = finite_value(finite_value(v) - c);
        double t = finite_value(s + y);
        c = finite_value((t - s) - y);
        s = t;
    }
    return s;
}

double neumaier_sum(const Vector& values) {
    double s = 0.0, c = 0.0;
    for (double v : values) {
        v = finite_value(v);
        double t = finite_value(s + v);
        double error = std::abs(s) >= std::abs(v)
            ? finite_value((s - t) + v) : finite_value((v - t) + s);
        c = finite_value(c + error);
        s = t;
    }
    return finite_value(s + c);
}

double pairwise_range(const Vector& values, std::size_t lo, std::size_t hi) {
    if (lo == hi) return 0.0;
    if (hi - lo == 1) return finite_value(values[lo]);
    std::size_t mid = lo + (hi - lo) / 2;
    return finite_value(pairwise_range(values, lo, mid) + pairwise_range(values, mid, hi));
}
double pairwise_sum(const Vector& values) { return pairwise_range(values, 0, values.size()); }

Vector solve_pivoted(const Matrix& a, const Vector& b, double pivot_tol = 1e-14) {
    std::size_t n = b.size();
    if (n == 0 || a.size() != n) throw std::invalid_argument("expected a nonempty square system");
    if (!std::isfinite(pivot_tol) || !(pivot_tol > 0.0 && pivot_tol < 1.0))
        throw std::invalid_argument("invalid pivot tolerance");
    Matrix rows(n, Vector(n + 1));
    double scale = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        if (a[i].size() != n) throw std::invalid_argument("expected a square matrix");
        for (std::size_t j = 0; j < n; ++j) {
            rows[i][j] = finite_value(a[i][j]);
            scale = std::max(scale, std::abs(rows[i][j]));
        }
        rows[i][n] = finite_value(b[i]);
    }
    if (scale == 0.0) throw std::invalid_argument("zero matrix");
    for (std::size_t j = 0; j < n; ++j) {
        std::size_t p = j;
        for (std::size_t i = j + 1; i < n; ++i)
            if (std::abs(rows[i][j]) > std::abs(rows[p][j])) p = i;
        if (std::abs(rows[p][j]) / scale <= pivot_tol)
            throw std::invalid_argument("pivot too small for this scale");
        std::swap(rows[j], rows[p]);
        for (std::size_t i = j + 1; i < n; ++i) {
            double factor = finite_value(rows[i][j] / rows[j][j]);
            rows[i][j] = 0.0;
            for (std::size_t k = j + 1; k <= n; ++k)
                rows[i][k] = finite_value(rows[i][k] - factor * rows[j][k]);
        }
    }
    Vector x(n, 0.0);
    for (std::size_t ii = n; ii > 0; --ii) {
        std::size_t i = ii - 1;
        Vector terms;
        for (std::size_t j = i + 1; j < n; ++j) terms.push_back(finite_value(rows[i][j] * x[j]));
        double rhs = finite_value(rows[i][n] - neumaier_sum(terms));
        x[i] = finite_value(rhs / rows[i][i]);
    }
    return x;
}

double residual_inf(const Matrix& a, const Vector& b, const Vector& x) {
    std::size_t n = b.size();
    if (n == 0 || a.size() != n || x.size() != n) throw std::invalid_argument("incompatible residual dimensions");
    double result = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        if (a[i].size() != n) throw std::invalid_argument("incompatible residual dimensions");
        Vector terms = {finite_value(b[i])};
        for (std::size_t j = 0; j < n; ++j)
            terms.push_back(finite_value(-finite_value(a[i][j]) * finite_value(x[j])));
        result = std::max(result, std::abs(neumaier_sum(terms)));
    }
    return result;
}

int main() {
    if (!std::numeric_limits<double>::is_iec559 || std::numeric_limits<double>::digits != 53)
        throw std::runtime_error("this example requires binary64");
    double m = 9007199254740992.0;
    Vector values = {m, 1.0, 1.0, -m}, cancellation = {1e16, 1.0, -1e16};
    Vector results = {naive_sum(values), kahan_sum(values), neumaier_sum(values),
                      pairwise_sum(values), kahan_sum(cancellation), neumaier_sum(cancellation)};
    std::vector<std::string> labels = {"trace naive", "trace kahan", "trace neumaier", "trace pairwise",
                                        "cancellation kahan", "cancellation neumaier"};
    std::cout << std::fixed << std::setprecision(0);
    for (std::size_t i = 0; i < results.size(); ++i) std::cout << labels[i] << " = " << results[i] << '\n';
    Vector expected = {0., 2., 2., 1., 0., 1.};
    std::cout << "summation checks = " << (results == expected ? "PASS" : "FAIL") << '\n';
    Matrix a = { {1e-20, 1.0}, {1.0, 1.0} };
    Vector b = {1.0, 2.0}, x = solve_pivoted(a, b);
    std::cout << std::setprecision(6) << "pivoted x = [" << x[0] << ", " << x[1] << "]\n";
    std::cout << "residual check = " << (residual_inf(a, b, x) <= 1e-12 ? "PASS" : "FAIL") << '\n';
    return 0;
}

import java.util.Arrays;
import java.util.Locale;

public class Main {
    static double finite(double value) {
        if (!Double.isFinite(value)) throw new IllegalArgumentException("nonfinite input or intermediate result");
        return value;
    }

    static double naiveSum(double[] values) {
        double s = 0.0;
        for (double value : values) s = finite(s + finite(value));
        return s;
    }

    static double kahanSum(double[] values) {
        double s = 0.0, c = 0.0;
        for (double value : values) {
            double y = finite(finite(value) - c);
            double t = finite(s + y);
            c = finite((t - s) - y);
            s = t;
        }
        return s;
    }

    static double neumaierSum(double[] values) {
        double s = 0.0, c = 0.0;
        for (double value : values) {
            value = finite(value);
            double t = finite(s + value);
            double error = Math.abs(s) >= Math.abs(value)
                ? finite((s - t) + value) : finite((value - t) + s);
            c = finite(c + error);
            s = t;
        }
        return finite(s + c);
    }

    static double pairwiseRange(double[] values, int lo, int hi) {
        if (lo == hi) return 0.0;
        if (hi - lo == 1) return finite(values[lo]);
        int mid = lo + (hi - lo) / 2;
        return finite(pairwiseRange(values, lo, mid) + pairwiseRange(values, mid, hi));
    }
    static double pairwiseSum(double[] values) { return pairwiseRange(values, 0, values.length); }

    static double[] solvePivoted(double[][] a, double[] b, double pivotTol) {
        int n = b.length;
        if (n == 0 || a.length != n) throw new IllegalArgumentException("expected a nonempty square system");
        if (!Double.isFinite(pivotTol) || !(pivotTol > 0.0 && pivotTol < 1.0))
            throw new IllegalArgumentException("invalid pivot tolerance");
        double[][] rows = new double[n][n + 1];
        double scale = 0.0;
        for (int i = 0; i < n; ++i) {
            if (a[i].length != n) throw new IllegalArgumentException("expected a square matrix");
            for (int j = 0; j < n; ++j) {
                rows[i][j] = finite(a[i][j]);
                scale = Math.max(scale, Math.abs(rows[i][j]));
            }
            rows[i][n] = finite(b[i]);
        }
        if (scale == 0.0) throw new IllegalArgumentException("zero matrix");
        for (int j = 0; j < n; ++j) {
            int p = j;
            for (int i = j + 1; i < n; ++i)
                if (Math.abs(rows[i][j]) > Math.abs(rows[p][j])) p = i;
            if (Math.abs(rows[p][j]) / scale <= pivotTol)
                throw new IllegalArgumentException("pivot too small for this scale");
            double[] swap = rows[j]; rows[j] = rows[p]; rows[p] = swap;
            for (int i = j + 1; i < n; ++i) {
                double factor = finite(rows[i][j] / rows[j][j]);
                rows[i][j] = 0.0;
                for (int k = j + 1; k <= n; ++k)
                    rows[i][k] = finite(rows[i][k] - factor * rows[j][k]);
            }
        }
        double[] x = new double[n];
        for (int i = n - 1; i >= 0; --i) {
            double[] terms = new double[n - i - 1];
            for (int j = i + 1; j < n; ++j) terms[j - i - 1] = finite(rows[i][j] * x[j]);
            double rhs = finite(rows[i][n] - neumaierSum(terms));
            x[i] = finite(rhs / rows[i][i]);
        }
        return x;
    }

    static double residualInf(double[][] a, double[] b, double[] x) {
        int n = b.length;
        if (n == 0 || a.length != n || x.length != n) throw new IllegalArgumentException("incompatible residual dimensions");
        double result = 0.0;
        for (int i = 0; i < n; ++i) {
            if (a[i].length != n) throw new IllegalArgumentException("incompatible residual dimensions");
            double[] terms = new double[n + 1]; terms[0] = finite(b[i]);
            for (int j = 0; j < n; ++j) terms[j + 1] = finite(-finite(a[i][j]) * finite(x[j]));
            result = Math.max(result, Math.abs(neumaierSum(terms)));
        }
        return result;
    }

    public static void main(String[] args) {
        double m = 9007199254740992.0;
        double[] values = {m, 1.0, 1.0, -m}, cancellation = {1e16, 1.0, -1e16};
        double[] results = {naiveSum(values), kahanSum(values), neumaierSum(values),
                            pairwiseSum(values), kahanSum(cancellation), neumaierSum(cancellation)};
        String[] labels = {"trace naive", "trace kahan", "trace neumaier", "trace pairwise",
                            "cancellation kahan", "cancellation neumaier"};
        for (int i = 0; i < results.length; ++i) System.out.printf(Locale.ROOT, "%s = %.0f%n", labels[i], results[i]);
        double[] expected = {0., 2., 2., 1., 0., 1.};
        System.out.println("summation checks = " + (Arrays.equals(results, expected) ? "PASS" : "FAIL"));
        double[][] a = { {1e-20, 1.0}, {1.0, 1.0} };
        double[] b = {1.0, 2.0}, x = solvePivoted(a, b, 1e-14);
        System.out.printf(Locale.ROOT, "pivoted x = [%.6f, %.6f]%n", x[0], x[1]);
        System.out.println("residual check = " + (residualInf(a, b, x) <= 1e-12 ? "PASS" : "FAIL"));
    }
}

All three programs print exactly these nine lines:

• trace naive = 0
• trace kahan = 2
• trace neumaier = 2
• trace pairwise = 1
• cancellation kahan = 0
• cancellation neumaier = 1
• summation checks = PASS
• pivoted x = [1.000000, 1.000000]
• residual check = PASS

The summation check verifies the expected successes and the deliberately exposed failures; it does not certify either method for all inputs. Empty sums return $0$. Nonfinite data and overflow are rejected. Empty or malformed systems, invalid pivot tolerances, zero matrices, and pivots below the guard are rejected. Inputs remain unchanged on both successful and rejected solves.

Complexity, Correctness, and Invariants

For a successful scan over $n$ stored values, naive, Kahan, and Neumaier summation all take $\Theta(n)$ time in the best, average, and worst case, with $\Theta(1)$ auxiliary state. There is no early termination based on the signs or sizes of ordinary finite inputs. Pairwise summation also takes $\Theta(n)$ work, with $\Theta(\log n)$ recursive stack space for $n\ge2$ in these implementations. Empty inputs take constant time and space. Rejecting an invalid input early can shorten a failed run.

For context, let $\gamma_k=ku/(1-ku)$ when $ku\lt1$, with $\gamma_0=0$. Under the usual rounding model with no overflow or problematic underflow, a naive nonempty sum satisfies

$$|\widehat S-S|\le\gamma_{n-1}\sum_{i=1}^{n}|x_i|.$$

The balanced tree replaces the long dependency chain by depth at most $d=\lceil\log_2 n\rceil$, giving the bound $\gamma_d\sum_i|x_i|$. A standard Kahan bound has leading behavior

$$|\widehat S-S|\le\bigl(2u+O(nu^2)\bigr)\sum_i|x_i|.$$

The last expression is an asymptotic error estimate, not a numerical certificate obtained by dropping the hidden higher-order term. Its small-leading-error interpretation assumes $nu$ is small and the required floating-point model holds. Dividing such bounds by $|S|$ introduces $\kappa_{\mathrm{sum}}$, so even a small accumulation coefficient may permit a large relative error on a strongly cancelling sum. The worked trace's exact Kahan answer is much better than the conservative bound predicts.

The naive summation bound follows by expanding the computation into a sum of terms multiplied by products of rounding factors. Each input passes through at most $n-1$ rounded additions after its exact addition to the initial zero; the accumulated factor differs from $1$ by at most $\gamma_{n-1}$. In a balanced tree each input instead passes through at most $d$ additions, giving the corresponding bound. These are worst-case guarantees, not predictions of typical signs of error.

For Kahan, it is tempting to claim that $s-c$ always equals the exact prefix sum. That equality is false in general: forming $y$ and computing the correction can also round. A more careful accounting makes the mechanism explicit. Define exact real error quantities

$$\rho_i=y_i-(x_i-c_{i-1}),\qquad \eta_i=t_i-(s_{i-1}+y_i),\qquad \zeta_i=c_i-\eta_i.$$

Because $s_i=t_i$, substitution gives the exact identity

$$s_i-c_i=(s_{i-1}-c_{i-1})+x_i+\rho_i-\zeta_i.$$

Thus the main addition error $\eta_i$ is cancelled in this accounting, leaving the errors from forming the adjusted term and recovering the correction. Establishing the stated bound requires controlling those remaining errors and the final correction as well. This explains both the improvement and why treating compensation as an exact extra precision register is wrong.

For the linear solve, the exact-arithmetic elimination invariant is that swapping rows and subtracting a multiple of one equation from another preserve the solution set. Floating-point elimination only approximately preserves it. Partial pivoting makes each active-column multiplier have magnitude at most $1$, while the growth of other entries still matters; there is no unconditional guarantee that every pivoted factorization has tiny error.

A scale-aware residual diagnostic for $\widehat x$ is

$$\eta_{\mathrm{norm}}=\frac{\|b-A\widehat x\|_\infty}{\|A\|_\infty\|\widehat x\|_\infty+\|b\|_\infty},$$

where $\|A\|_\infty=\max_i\sum_j|a_{ij}|$. Treat an exactly zero numerator and denominator as zero; large scales may require careful evaluation to avoid overflow. This quantity measures normwise backward error when both $A$ and $b$ may be perturbed. It does not measure forward error directly. When $A$ is nonsingular and $b\ne0$, the residual relation gives

$$\frac{\|\widehat x-x\|_\infty}{\|x\|_\infty}\le\kappa_\infty(A)\frac{\|b-A\widehat x\|_\infty}{\|b\|_\infty},\qquad\kappa_\infty(A)=\|A\|_\infty\|A^{-1}\|_\infty.$$

Consequently, a large condition number can turn a very small relative residual into an appreciable solution error. The demo prints only its simpler absolute residual test, and does not estimate a condition number or certify this bound numerically.

Variants, Pitfalls, and Optimizations

Compensation improves a particular source of error. Choose the response that matches the actual failure:

• Reformulate a cancellation-prone expression. For positive $x$, replace $\sqrt{x+1}-\sqrt{x}$ by $1/(\sqrt{x+1}+\sqrt{x})$. At stored $x=10^{16}$, the first expression returns $0$ because $x+1$ rounds to $x$; the second yields about $5\times10^{-9}$. This helps evaluate the function on the stored $x$ without relying on the small difference of two rounded square roots.
• Use specialized elementary functions. For small $z$, log1p(z) and expm1(z) avoid first forming a rounded $1+z$ or subtracting $1$ from a rounded exponential. Their value is in the formulation, not merely in using a wider type.
• Choose a reduction suited to the data and execution model. Kahan suits streaming input; Neumaier handles the displayed late-cancellation case better. Balanced pairwise reduction fits parallel trees. Combining only the main totals from compensated chunks discards their corrections; use a specified merge algorithm or exact accumulator when the guarantee requires it.
• Separate summation error from product error. A compensated dot product that first forms $\operatorname{fl}(a_i b_i)$ cannot restore bits already lost in those products. Fused multiply-add can help product-error recovery when used deliberately in a designed algorithm, but contracting arbitrary expressions changes the arithmetic being analyzed.
• Use a suitable representation. Integer minor units or decimal arithmetic may fit fixed-scale financial quantities; arbitrary precision or exact rational arithmetic can provide references. Increasing binary precision postpones failures but does not remove conditioning or justify an unstable formula.
• Check linear algebra in the original problem. Pivoting, scaling, a residual evaluated with stronger arithmetic, and sometimes iterative refinement address different issues. Compensating a residual cannot undo a bad elimination step, and partial pivoting does not fix ill-conditioning.

More implementation traps are easy to miss. Updating $s$ before computing $c$ destroys Kahan's use of the old total. Optimizing $(t-s)-y$ to zero by treating real algebra as floating-point algebra removes compensation. A universal comparison such as abs(a-b) < 1e-9 ignores scale; a practical test often combines an application-specific absolute tolerance with a relative tolerance, for example $|a-b|\le\mathrm{atol}+\mathrm{rtol}\max(|a|,|b|)$. Neither a tolerance nor a small residual can replace a trustworthy reference when validating numerical accuracy.

Changing the order of inputs or a parallel reduction tree can change the answer. A deterministic result may still be inaccurate, and a more accurate result may differ bitwise from an earlier version. Specify whether the requirement is an error bound, a correctly rounded answer, or repeatability under a particular order and environment.

Empty input, signed zero, mixed signs, huge exponent gaps, subnormals, nonfinite data, and sums that overflow despite a finite exact final total deserve deliberate tests. These implementations reject nonfinite intermediate sums rather than rearranging them to avoid overflow. They accept finite subnormal values, but do not certify the normal-range error model for computations involving underflow. Signed-zero details are not part of their correctness contract.

Practice Problems and Real-World Use

For a meaningful summation test, construct an exact reference from the actual stored values, not from their short printed decimal spellings. In Python, convert each float with as_integer_ratio() or Fraction.from_float, add the resulting rationals exactly, and compare the computed binary64 result as another exact rational. Use absolute error when the true sum is zero and report relative error only when it is defined. A high-precision decimal reference must likewise distinguish decimal intentions from exact conversions of the stored binary values.

Test the expected failure cases as well as the successes. Include all permutations of small cancelling multisets, the worked trace, the Neumaier example, empty input, singletons, positive-only lists, and reproducible random values across a wide exponent range. Record each method's actual error; do not assert that Kahan always beats naive summation or that pairwise summation is always exact. For matrices, compare with an independent solve, verify the inputs were preserved, and evaluate the residual in stronger arithmetic.

For further reading, Goldberg's floating-point arithmetic tutorial discusses rounding, cancellation, and compensated summation. The Python mathematical-functions documentation describes math.fsum, which tracks multiple partial sums, and log1p and expm1. A library routine's documented accuracy and exceptional-input behavior should guide production use; math.fsum is a useful comparison here, while exact rational arithmetic provides the independent stored-input reference.