# Proof that (x-y)(x+y) is more accurate than x²-y²

I was carrying on my reading of What Every Computer Scientist Should Know About Floating-Point Arithmetic but got stuck on the proof of Theorem 2 (page 34).

At some point it says: \begin{align} (x \otimes x) \ominus (y \otimes y) & = \left[x^2(1 + \delta_1) - y^2(1 + \delta_2)\right](1 + \delta_3) \\ & = \left[(x^2 - y^2)(1 + \delta_1) + (\delta_1 - \delta_2)y^2\right](1 + \delta_3) \\ \end{align} I'm ok with the rewriting, but I don't understand the argument that:

When $$x$$ and $$y$$ are nearby, the error term $$(\delta_1 - \delta_2)y^2$$ can be as large as the result $$x^2 - y^2$$.

This doesn't make much sens to me right now. I understand that if both quantities are close to each other, then the relative error is close to $$1$$. But not why they should be close to each other.

Suppose $$y = x + \Delta$$. Then $$(x^2 - y^2) = x^2 - (x^2 + 2\Delta x + \Delta^2) = - (2\Delta x + \Delta^2)$$ with leading term on the order of $$2 \Delta x$$. Multiply by $$1 + \delta_1$$ and that's still the leading term.

Compare to $$(\delta_1 - \delta_2) y^2$$ with leading term $$(\delta_1 - \delta_2) x^2$$.

I'm not 100% sure of everything, but here are some elements that are implicit in the proof. Note that I reuse the notations and common knowledge of the paper, therefore this answer isn't self-contained.

First, let's assume that $$x \gt y \ge 0$$. It's a safe assumption to make because the sign of the variables doesn't matter and swapping $$x$$ and $$y$$ only change the sign of the result.

Let's write the relative error for the formula $$(x \otimes x) \ominus (y \otimes y)$$.

\begin{align} relerr & = \frac{\left|(x \otimes x) \ominus (y \otimes y) - (x^2 - y^2)\right|}{\left|x^2 - y^2\right|} \\ & = \left|\frac{(x \otimes x) \ominus (y \otimes y) - (x^2 - y^2)}{x^2 - y^2}\right| \\ & = \left|\frac{(x \otimes x) \ominus (y \otimes y)}{x^2 - y^2} - 1\right| \\ & = \left|\frac{\left[(x^2 - y^2)(1 + \delta_1) + (\delta_1 - \delta_2)y^2\right](1 + \delta_3)}{x^2 - y^2} - 1\right| \\ & = \left|\frac{(\delta_1 - \delta_2)y^2}{x^2 - y^2}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ \end{align}

Since we want to bound this expression by some amount of the machine-$$\epsilon$$, we need to find when the first term is the largest. Which happen to be when $$y$$ is large, making the numerator large, and when $$x$$ is close to $$y$$, making the denominator small. We could probably stop there (as the paper does), but this doesn't provide an upper bound for the relative error.

Let's assume that $$x = (1 + \delta_x) y$$ with $$0 < \delta_x \ll 1$$.

\begin{align} relerr & = \left|\frac{(\delta_1 - \delta_2)y^2}{(1 + \delta_x)^2 y^2 - y^2}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ & = \left|\frac{\delta_1 - \delta_2}{(1 + \delta_x)^2 - 1}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ & = \left|\frac{\delta_1 - \delta_2}{2 \delta_x + \delta_x^2}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ \end{align}

In order to bound this expression, we will need a lower bound on $$\delta_x$$. It needs to be greater than $$0$$ so that $$x \neq y$$ and the equations make sens. The smallest non-zero value it can take is $$\delta_x \ge \beta^{-p+1} = 2 \epsilon$$.

Moreover we have: $$-\epsilon \le \delta_1 \le \epsilon$$ and $$-\epsilon \le \delta_2 \le \epsilon$$.

Therefore, we have: \begin{align} \delta_1 - \delta_2 & \le 2 \epsilon \\ 2 \delta_x + \delta_x^2 & \ge 4 \epsilon + 4 \epsilon^2 \\ \frac{\delta_1 - \delta_2}{2 \delta_x + \delta_x^2} & \le \frac{2 \epsilon}{4 \epsilon + 4 \epsilon^2} & = \frac{1}{2 + \epsilon} & \lt \frac{1}{2} \\ \end{align}

Hence: \begin{align} relerr & \le \left|\frac{\delta_1 - \delta_2}{2 \delta_x + \delta_x^2}(1 + \delta_3)\right| + \left|(1 + \delta_1)(1 + \delta_3) - 1\right| \\ & \lt \frac{1 + \delta_3}{2} + |\delta_1 + \delta_3 + \delta_1 \delta_3| \\ & \lt \frac{1 + 2 \epsilon}{2} + \epsilon + 2 \epsilon + 2 \epsilon^2 \\ & \lt \frac{1}{2} + 4 \epsilon + 2 \epsilon^2 \\ & \lt \frac{1}{2} + 5 \epsilon \end{align}

Meaning that we could get an answer between 50% and 150% of the true result. I was however unable to find an actual example by exhaustive search in the high numbers. The best I could find (or worse, depending how you look at it ^^) was: $$x = 1.340780783004664349 \times 10^{154}$$ $$y = 1.340780783004664051 \times 10^{154}$$ using 64 bit floats (53 bits mantissa). Which has a relative error of roughly $$0.25$$.

Given that the binary mantissa of these numbers looks pretty regular, there might be something special about them and my upper bound is not correct.

My guess is that when $$x$$ is the float right after $$y$$, then maybe the rounding error of $$x^2$$ is similar to that of $$y^2$$. Which would make $$\delta_1 - \delta_2 \le \epsilon$$ instead of $$2 \epsilon$$. If anyone can prove it, feel free to edit this answer.