Catastrophic cancellation (source code)

= Catastrophic cancellation
{title2=$\sqrt{1+x}-1=\dfrac{x}{\sqrt{1+x}+1}$}
{wiki}

<Catastrophic cancellation> occurs when <subtraction> of nearly equal approximations removes their matching leading digits and exposes earlier <roundoff error>. The <subtraction> can itself be exact on the stored operands: the information may already have been lost while obtaining them. Goldberg distinguishes this from cancellation of accurately known quantities, which can be harmless.https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html[Goldberg, 1991, §Cancellation]{ref}

For $x>0$, consider the single <real-valued function>
$$
f(x)=\sqrt{1+x}-1=\frac{x}{\sqrt{1+x}+1}.
$$
The <rationalization of a square-root difference> proves equality in exact <arithmetic>. The two evaluations use the same input $x$, but the direct version first forms a value near one and then subtracts one. The rearranged version keeps the small input in its <numerator> and divides by a <denominator> near two.

\Image[/floating-point-cancellation.png]
{title=Direct and rationalized evaluation of a small square-root difference against a high-precision reference for the same binary64 inputs}

For binary64 `x = 1e-16`, the rounded <sum> `1.0 + x` is one, so `sqrt(1.0 + x) - 1.0` is zero. The rationalized evaluation returns approximately $5\times10^{-17}$. A 100-digit Decimal calculation using the exact stored input gives relative errors of one and about $2.50\times10^{-17}$ respectively. At `x = 1e-12`, the direct result is about $5.0004445029\times10^{-13}$, with relative error $8.89\times10^{-5}$; the rearranged result is about $4.99999999999875\times10^{-13}$, with relative error $1.21\times10^{-17}$.

The plot uses `Decimal.from_float(x)` to preserve the binary input exactly, computes $x/(\sqrt{1+x}+1)$ at 100 decimal digits, and verifies the reference at 160 digits. Exact conversion of a stored binary value is documented by Python.https://docs.python.org/3.14/library/decimal.html#decimal.Decimal.from_float[Python 3.14 documentation, Decimal.from_float]{ref} Each error is $|\widehat f-f|/|f|$, with the computed result converted exactly to Decimal before comparison. Hence the curves measure evaluation error for one common input, rather than a mixture of input conversion error and evaluation error. The plotted interval $10^{-16}\leq x\leq10^{-1}$ is far above the <subnormal number> range.

The original problem is well conditioned with respect to relative changes in $x$. Using the <derivative> $f'(x)=1/(2\sqrt{1+x})$ and <rationalization of a square-root difference>, its relative <condition number> is
$$
\kappa_f(x)=\left|\frac{xf'(x)}{f(x)}\right|=\frac{\sqrt{1+x}+1}{2\sqrt{1+x}},\qquad \frac12<\kappa_f(x)<1\quad(x>0).
$$
It tends to one as $x$ tends to zero. The direct evaluation instead introduces the intermediate difference $s-1$ with $s=\sqrt{1+x}$; an error of order <unit roundoff> in $s$ is large relative to $f(x)\approx x/2$. This is a failure of the chosen <algorithm>'s accuracy despite a small <condition number> for the same <real-valued function>. The rationalized evaluation avoids that intermediate <catastrophic cancellation>, retaining relative accuracy near the <unit roundoff> scale throughout the plotted interval. Recognizing this structure can improve a calculation far more than merely printing extra digits.