Three floating-point failures in summation, cancellation and RANDU

Each section computes something in double precision and compares it with an exact answer. The references are in posts/refs/kahan-summation-horner-noise-and-randus-planes.py (output beside it): Python's math.fsum, which returns the correctly rounded sum of its arguments and which the reference checks against exact rational arithmetic (fractions) up to $N=10^4$; the exact value of $(x-1{)}^7$ in rationals; and integer arithmetic for the generators.

A double has a 53-bit significand, so a sum, product or quotient of two doubles is the exact result times $(1+\delta )$ with $|\delta |\le 2^{-53}$. $\epsilon =2^{-52}$ is the gap between 1 and the next double, and one ulp (unit in the last place) of a number $s$ is the gap between $s$ and the next double above it in magnitude, about $\epsilon |s|$.

Summing $1/k$

Adding term by term, each addition rounds the running sum, and the rounding errors pile up. Once the sum is near 14, a term like $10^{-6}$ keeps only the bits that line up with the sum's last bit, about 30 of its 53. Kahan's method (1965) keeps a second variable $c$ holding what the last addition rounded away, and subtracts it from the next term before adding:

In exact arithmetic $c$ would be 0. In floating point, $(t-s)$ is the part of $y$ that made it into $t$, and subtracting $y$ leaves the part that did not, with its sign flipped. The error of the result is at most about $2\epsilon |s|$, independent of the number of terms, as long as $N\epsilon$ is small.

The cell sums $1/k$ (each term the double nearest $1/k$) up to $N=10^6$ both ways and records both sums at 31 checkpoints spaced evenly in $\log N$.

$(x-1{)}^7$ near its root

Expanded, $(x-1{)}^7=x^7-7x^6+21x^5-35x^4+35x^3-21x^2+7x-1$. Horner's rule evaluates it with seven multiplications and seven additions, $((((((x-7)x+21)x-35)x+35)x-21)x+7)x-1$. Near $x=1$ the partial results are 1, $-6$, 15, $-20$, 15, $-6$, 1, while the answer is tiny ($10^{-14}$ at $x=1.01$). Each partial result is rounded, by up to half an ulp of a number near 20, about $2\times 10^{-15}$, and those errors are as large as the answer. The rounding bound for Horner's rule (Higham, Accuracy and Stability of Numerical Algorithms, section 5.1) is ${\gamma }_{14}{\sum }_i|c_i||x{|}^i$ with ${\gamma }_{14}=14u/(1-14u)$ and the unit roundoff $u=2^{-53}$; here it is about $2\times 10^{-13}$, and the measured errors reach $2.6\%$ of it.

The factored form has no such problem. For $x$ between $1/2$ and 2, $x-1$ is computed exactly (Sterbenz's lemma: the difference of two doubles within a factor of 2 of each other is a double), and then six multiplications each add a relative error of at most $2^{-53}$.

At $x=1.2$, away from the root, Horner's value has a relative error of $6\times 10^{-11}$. Near the root it has no correct digits.

RANDU's planes

RANDU was IBM's random number generator in the 1960s: $x_{k+1}=65539x_k\mathrm{mod}{2}^{31}$, with draws $u_k=x_k/2^{31}$. Its multiplier is $2^{16}+3$, and squaring it gives $2^{32}+6\cdot 2^{16}+9$. The first term vanishes modulo $2^{31}$, and $6\cdot 2^{16}=6\cdot 65539-18$, so $65539^2\equiv 6\cdot 65539-9$, and every three consecutive values satisfy

So $9u_k-6u_{k+1}+u_{k+2}$ is always a whole number. With each $u$ in $[0,1)$ that number lies between $-5$ and 9, which puts every point $(u_k,u_{k+1},u_{k+2})$ on one of the 15 parallel planes $9u-6v+w=j$. A simulation that uses triples of draws as points in a cube samples 15 planes instead of the cube.

The products here are below $2^{48}$, so double-precision arithmetic computes them, and the remainders, exactly. The cell also runs MINSTD, the Park-Miller generator $x_{k+1}=16807x_k\mathrm{mod}(2^{31}-1)$, for which the same combination is never a whole number in 3000 triples.