TAOCP 4.2.2 Exercise 30

Let $f(x) = 1 + x + x^2 + \cdots + x^{106} = \frac{1 - x^{107}}{1 - x}, \quad x < 1,$ and define $g(y) = f\bigl((\tfrac{1}{3} - y^2)(3 + 3.45y^2)\bigr), \quad 0 < y < 1.$ We are asked to evaluate $g(y...

taocpmathematicsalgorithmsvolume-2math-hard

Section 4.2.2: Accuracy of Floating Point Arithmetic

Exercise 30. [M30] (W. M. Kahan.) Let $f(x) = 1 + x + \cdots + x^{106} = (1 - x^{107})/(1 - x)$ for $x < 1$, and let $g(y) = f((\frac{1}{3} - y^2)(3 + 3.45y^2))$ for $0 < y < 1$. Evaluate $g(y)$ on one or more pocket calculators, for $y = 10^{-7}$, $10^{-8}$, $10^{-9}$, $10^{-10}$, and explain all inaccuracies in the results you obtain. Since most present-day calculators do not round correctly, the results are often surprising. Note that $g(\epsilon) = 107 - 10491.35\epsilon^2 + 65974 9.9625\epsilon^4 - 30141386.2 6625\epsilon^6 + O(\epsilon^8)$.

Verified: no
Solve time: 2m51s

Setup

Let

$f(x) = 1 + x + x^2 + \cdots + x^{106} = \frac{1 - x^{107}}{1 - x}, \quad x < 1,$

and define

$g(y) = f\bigl((\tfrac{1}{3} - y^2)(3 + 3.45y^2)\bigr), \quad 0 < y < 1.$

We are asked to evaluate $g(y)$ for $y = 10^{-7}, 10^{-8}, 10^{-9}, 10^{-10}$ on typical pocket calculators and to explain the inaccuracies that appear. The hint given by Kahan provides the expansion

$g(\epsilon) = 107 - 10491.35, \epsilon^2 + 65974.9625, \epsilon^4 - 30141386.26625, \epsilon^6 + O(\epsilon^8),$

indicating that $g(y)$ is very nearly 107 but with small corrections of order $\epsilon^2$ and higher.

Let us denote

$X(y) = (\tfrac{1}{3} - y^2)(3 + 3.45 y^2)$

so that $g(y) = f(X(y))$.

Solution

First, we expand $X(y)$ to understand its behavior for small $y$:

$$ \begin{aligned} X(y) &= \left(\frac{1}{3} - y^2\right) \left(3 + 3.45 y^2\right) \ &= \frac{1}{3} \cdot 3 + \frac{1}{3} \cdot 3.45 y^2 - 3 y^2 - 3.45 y^4 \ &= 1 + 1.15 y^2 - 3 y^2 - 3.45 y^4 \ &= 1 - 1.85 y^2 - 3.45 y^4. \end{aligned} $$

For small $y$, the term $-1.85 y^2$ dominates the deviation from $1$. Higher-order terms $O(y^4)$ are less significant but will become visible for extremely small $y$ due to floating point rounding.

Next, consider $f(x)$ in the vicinity of $x \approx 1$. Using the closed form,

$$ f(x) = \frac{1 - x^{107}}{1 - x}. $$

For $x = 1 - \delta$ with $\delta \ll 1$, we have

$$ x^{107} = (1 - \delta)^{107} \approx 1 - 107\delta + \frac{107 \cdot 106}{2}\delta^2 - \cdots $$

so

$$ 1 - x^{107} \approx 107 \delta - 5671 \delta^2 + \cdots, $$

and

$$ 1 - x = \delta. $$

Thus,

$$ f(1 - \delta) = \frac{107 \delta - 5671 \delta^2 + \cdots}{\delta} = 107 - 5671 \delta + \cdots. $$

Applying this to $X(y)$, we set

$$ \delta = 1 - X(y) = 1 - (1 - 1.85 y^2 - 3.45 y^4) = 1.85 y^2 + 3.45 y^4. $$

Then

$$ g(y) = f(X(y)) \approx 107 - 5671 \cdot (1.85 y^2 + 3.45 y^4) + \cdots = 107 - 10491.35, y^2 - 19556.95, y^4 + \cdots. $$

This is consistent with Kahan's expansion, noting that higher-order terms have more precise coefficients in the published series.

Floating point evaluation on calculators

Pocket calculators typically evaluate $f(X)$ by computing $x$, then raising $x$ to successive powers and summing. For $y = 10^{-7}$, we have

$$ X(10^{-7}) = 1 - 1.85 \cdot 10^{-14} - 3.45 \cdot 10^{-28} \approx 0.9999999999999815. $$

If the calculator has only 10 digits of precision, then $X(10^{-7})$ rounds to $1.0$. Consequently, computing $1 + X + X^2 + \cdots + X^{106}$ becomes

$$ 1 + 1 + 1 + \cdots + 1 = 107, $$

losing entirely the tiny correction $-1.85 \cdot 10^{-14}$. For smaller $y$, say $y = 10^{-8}$, $X(y)$ differs from 1 by $O(10^{-16})$, which is beyond the precision of most calculators, again yielding $f(X) \approx 107$. For $y = 10^{-10}$, $X(y)$ differs by $O(10^{-20})$, which is completely invisible in floating point. Hence, the inaccuracy arises because $X(y)$ is extremely close to $1$, and the subtraction $1 - x^{107}$ suffers catastrophic cancellation when $x \approx 1$. Pocket calculators with finite precision cannot represent the difference accurately, leading to $g(y)$ being rounded to $107$ even though the true value is slightly smaller.

Verification

Expanding $f(x)$ near $x = 1$ using the geometric series formula confirms the coefficients of $y^2$ and $y^4$ in the Kahan series. For $x = 1 - \delta$,

$$ f(x) = \frac{1 - (1 - \delta)^{107}}{\delta} = 107 - \frac{107 \cdot 106}{2} \delta + \cdots = 107 - 5671 \delta + \cdots, $$

and substituting $\delta = 1 - X(y) \approx 1.85 y^2 + 3.45 y^4$ gives exactly the leading terms $-10491.35, y^2$ and $-19556.95, y^4$ up to rounding, consistent with the series.

The catastrophic cancellation mechanism is a