# how do you compute $\|c-a\| - \|b-a\|$ without catastrophic cancellation?

Given three points or vectors in the plane: \begin{align} \vec a &= (a_x,a_y) \\ \vec b &= (b_x,b_x) \\ \vec c &= (c_x,c_y) \end{align} How do you compute $\lVert \vec c - \vec a \rVert - \lVert \vec b - \vec a \rVert$, i.e. "how much farther is it from $\vec a$ to $\vec c$ than from $\vec a$ to $\vec b$"?

For definiteness, all coordinates and computatations are to be in double precision IEEE754 floating point arithmetic. The answer must be reasonably accurate even if $\vec a$ is very large compared to $\vec b$ and $\vec c$.

Note that the naive expression $$\sqrt{{(c_x-a_x)}^2+{(c_y-a_y)}^2} - \sqrt{{(b_x-a_x)}^2+{(b_y-a_y)}^2},$$ while mathematically correct, is unsuitable for this computation because it catastrophically cancels if $\vec a$ has much greater magnitude than $\vec b$ and $\vec c$.

For example, if: \begin{align} \vec a &= (-10^{20},-10^{20}) \\ \vec b &= (0,0) \\ \vec c &= (1,1) \end{align} then the answer is $\sqrt2$, but computing it naively will produce $0$ due to catastrophic cancellation.

• Conjugate, probably
– Pedro
Dec 30, 2016 at 11:30
• @PedroTamaroff Yes, probably, but it doesn't seem to be that simple. Dec 30, 2016 at 11:43
• I've never heard the phrase "catastrophic cancellation" before, but I'll be sure to work it into conversation the next chance I get.
– Paul
Dec 30, 2016 at 13:06

The binomial identity $x^2-y^2=(x+y)(x-y)$ can be extended to euclidean norms and scalar products, which gives $$\|c-a\|-\|b-a\|=\frac{\|c-a\|^2-\|b-a\|^2}{\|c-a\|+\|b-a\|} =\frac{\bigl\langle b+c-2a,\,c-b\bigr\rangle}{\|c-a\|+\|b-a\|}$$
• To belabor the second point in my previous comment: I finally got clarity when I realized that you actually used that binomial identity pattern (aka conjugation) twice: first the scalar form $x^2-y^2=(x+y)(x-y)$ with $x=||\vec c - \vec a||$ and $y=||\vec b - \vec a||$, and then the vector form $<\vec x,\vec x>-<\vec y,\vec y>=<\vec x+\vec y,\vec x-\vec y>$ with $\vec x=\vec c-\vec a$ and $\vec y=\vec b-\vec a$. This gives two ways out from your middle expression; choosing the wrong one (due to not realizing there was a choice!) kept leading me back to where I started, at first. Dec 31, 2016 at 4:49