Why is rewriting $x^2 -y^2$ as $(x+y)(x-y)$ a way to avoid catastrophic cancellation?
We are still doing $(x-y)$. Is it because the last operation in the second form is a multiplication?
[Math] How does rewriting $x^2 -y^2$ as $(x+y)(x-y)$ avoid catastrophic cancellation
catastrophic-cancellationfloating pointnumerical methods
Related Solutions
One way to rewrite the expression is to multiply by the conjugate:
$$\frac{1 - \sqrt{1 + 5x}}{x} \cdot \frac{1 + \sqrt{1 + 5x}}{1 + \sqrt{1 + 5x}} = \frac{1 - (1 + 5x)}{x(1 + \sqrt{1 + 5x})} = -\frac{5}{1 + \sqrt{1 + 5x}}$$
For this second expression, when $x \sim 0$ the denominator is $\sim 2$, and you get a number close to $-5/2$. In the original expression when $x \sim 0$ both the numerator and denominator are close to zero so you could get problems.
But it doesn't actually seem like you got any problems in your computation, so...?
When $\cos(x)$ is computed for a reasonably small $x$, a number like 0.99999999999994598649857etc is obtained; the IEEE double precision standard defines a way to round this number to roughly 16 decimals. Thus out comes something like 0.99999999999995. I didn't count the digits, but you get the idea (and anyway, rounding is done in base 2). This does not seem catastrophic. But now subtract 1 from both numbers, and you see that a huge relative error has been committed to that result. The division by $x^2$ roughly gives you this relative error, because when $x$ is small, then $\cos(x) \approx 1 - \frac12 x^2$ (according to Taylor series). So when $x$ is very small, $\cos(x)$ is rounded to 1 (no digits after the comma), which explains the plateau in your plot.
Follow up on OP's request: To formalize this, let $fl$, for "floating point arithmetics", be the function that rounds a real number to the nearest machine number ($=$ a number representable as double precision IEEE number, of which there are only finitely many). Although maybe not strictly accurately speaking (there may be optimization; there may be hidden extended precision computation; and the rounding is not always to the nearest machine number), the machine computes not $(1 - \cos(x)) / x^2$ but $$ fl( fl( 1 - {\color{red}{fl}}(\cos(x)) ) / fl(x^2) ) , $$ supposing we start with $x$ such that $x = fl(x)$. The issue originates with the application of the $fl$ that is marked in red, because $$ {\color{red}{fl}}(\cos(x)) = \cos(x) (1 + \delta) $$ where $|\delta| \leq \text{eps} \approx 10^{-16}$, and $\text{eps}$ is the smallest positive number such that $$ fl(1 + \text{eps}) \neq 1 . $$ Therefore, $$ 1 - {\color{red}{fl}}(\cos(x)) = 1 - \cos(x) + \delta \cos(x) \approx \tfrac12 x^2 + \delta $$ when $x$ is not too small (otherwise ${\color{red}{fl}}(\cos(x)) = 1$). When $x$ is of order $\sqrt{\text{eps}} \approx 10^{-8}$ (as in your picture), this result $\tfrac12 x^2$ is contaminated by an absolute error of $\delta$, or a relative error of $2 \delta / x^2$. So the final result $(1 - \cos(x)) / x^2$ is contaminated by an absolute error of $\delta / x^2 \approx 1$ when $x \approx \sqrt{\text{eps}}$, which, again, is what you see in the picture.
ps. There have been infamous consequences of finite precision arithmetics.
Best Answer
$x^2-y^2$ will be only as precise as $x^2$ is. If $x$ and $y$ are close to each other, the computational error might be larger than the result.
In the second case, less loss of significance occurs.
Consider e.g. $x = 9000.2$ and $y = 9000.1$, where the floating-point precision is 5 significant digits. $x^2-y^2$ is actually equal to 1800.03.
Let's try to compute it in $x^2-y^2$ way: $x^2 = 81003600.04 = 81003000$; $y^2 = 81001800.01 = 81001000$, so $x^2-y^2 = 2000$, which is quite far away from the precise answer.
So let's try to compute it in $(x+y)(x-y)$ way: $x+y = 18000.3 = 18000$; $x-y = 0.1$, so $(x+y)(x-y) = 1800$, which is significantly closer to the precise answer (and is actually as close as possible with 5 significant digits).