6

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.

babakks
  • 1,660
Don Hatch
  • 1,277

1 Answers1

8

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\|} $$

Don Hatch
  • 1,277
Lutz Lehmann
  • 131,652
  • 1
    This of course is the explicit demonstration of Pedro's hint in the comments. – J. M. ain't a mathematician Dec 30 '16 at 12:49
  • This solves my problem very nicely; thank you. A couple of things: (1) Is there a reason you negated everything wrt to my original formulation? The negation makes it a bit harder for me to see the geometric meaning of both the initial expression and the answer. (2) It took me a while to see how you arrived at the final equality from the binomial identity. In my own notes on how to do this, I'm adding two intermediate steps, so the whole thing becomes: ||c-a||-||b-a|| = (||c-a||^2-||b-a||^2)/(...) = (<c-a,c-a>-<b-a,b-a>)/(...) = <(c-a)+(b-a),(c-a)-(b-a)>/(...) = <b+c-2a,c-b>/(...) – Don Hatch Dec 30 '16 at 23:40
  • 1
    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. – Don Hatch Dec 31 '16 at 04:49
  • For the negation there is no reason, I just started without checking and then there was too much to change. But as you saw, the difference in the end is the negation of both factors in the scalar product. – Lutz Lehmann Dec 31 '16 at 07:39
  • That's beautiful. – Danra Oct 31 '17 at 18:08