Numerically stable algorithm for the angle between two vectors

[Mincheol_L_]

• Is there anyone who knew the numerically stable algorithm for the angle between two vectors? 
• - atan2 is more stable than acos?
• - acos is more stable than atan2, if the two vectors have normalized?

[andrew_4619]

Are these 2D or 3D vectors and do you want the acute angle? atan2 gives you the answers in the -PI to  PI range, acos is 0 to PI and asin is  -pi/2 to pi/2 so it depends what you want.

I don't think you will find any difference in the stability, atan2 is the most versatile and informative. Use real(8) if accuracy beyond 7SF is needed.

For 3D I normalise the vectors , do an asin of the dot product.

[Mincheol_L_]

Thanks app4619.
I want both acute angle and obtuse angle.
The algorithm have to be different for 2D and 3D vectors?
If so, please give me more detailed information.

[andrew_4619 · Accepted solution]

If the vectors were 2D (a1,a2) & (b1,b2)  the angle would be atan2(b2,b1)-atan2(a2,a1) the angle being positive as the x axis rotates towards the y axis. The only problem case is if either of the vectors is 0,0 which is indeterminate (x,0) and (0,y) are OK.

for 3D you have acos( dot_product(v1,v2)/(norm2(v1)*norm2(v2)))  where v1 and v2 are 3d vectors.This one is always going to give an acute angle and again null vectors will kill it (div zero).

[Anthony_Richards]

How about

eta=+ve number a litle bigger than the inverse of the biggest +ve real number representable by the machine
norm1=x1*x1+y1*y1+z1*z1
if(norm1.lt.eta) ifail=-1, return
norm2=(x2*x2+y2*y2+z2*z2 )*norm1
if(norm2.lt.eta) ifail=-1, return
angle=acos( (x1*x2+y1*y2+z1*z2)/norm2 )
ifail=0, return

?

[jimdempseyatthecove]

Try to scale the numbers such that you do not have overflow/underflow

HowBig1 = abs(x1) + abs(y1) + abs(z1)
if(HowBig1 .eq. 0.0) HowBig1 = 1.0
HowBig2 = abs(x2) + abs(y2) + abs(z2)
if(HowBig2 .eq. 0.0) HowBig2 = 1.0
x1fixed = x1 / HowBig1
y1fixed = y1 / HowBig1
z1fixed = z1 / HowBig1
... same for the 2's
norm1 = x1fixed**2 + y1fixed**2 + z1fixed**2
norm2 = (x2fixed**2 + y2fixed**2 + z2fixed**2)*norm1
if(norm1.lt.eta) ifail=-1, return
angle=acos( (x1fixed*x2fixed+y1*y2fixed+z1fixed*z2fixed)/norm2 )

Jim Dempsey

[NotThatItMatters]

Where are the square roots in the norms?

[jimdempseyatthecove]

That is for you to add. The snip of code provided by Anthony did not include them so I did not wish to clutter my example with additional statements. The intent of the example was to illustrate how to scale the input such that the sum of the squares would not overflow.

I will point you in the direction, you can scan the code for omissions.

... do the fixed code from earlier
norm1 = sqrt(x1fixed**2 + y1fixed**2 + z1fixed**2)
norm2 = sqrt(x2fixed**2 + y2fixed**2 + z2fixed**2)
norm12 = norm1 * norm2
if(norm12.lt.eta) ifail=-1, return
angle=acos( (x1fixed*x2fixed+y1*y2fixed+z1fixed*z2fixed)/norm12 )

The ..fixed numbers are on the order of 0.3, the norms are on the order of sqrt(2.7) or 1.65, the inner expression of acos is on the order of 1.

Jim Dempsey

[andrew_4619]

norm2 is an F2008 instrinsic function BTW chaps....

[jimdempseyatthecove]

The code snip is abstract. The variable names chosen were those of the referenced snip from the prior post. Changing the name to avoid conflict may have confused the intent.  Length1, Length2, and Length1xLength2 might have been more appropriate choice of variable names. Also Fortran permits you do define a variable with the same name as an intrinsic function and keywords (REAL :: DO). I don't own a pair of chaps.

Overflow and Underflow remediation is a technique I am trying to teach. However, while what I've shown will reduce the probability of overflow/underflow it will affect the precision of the result as it adds operations, each with the potential of introducing additional round off error.

Mincheol will have to take these issues under consideration when crafting the code to perform the desired functionality (or find a robust enough library function that performs in the manner desired).

BTW your #4 suggestion is likely the way to go... unless norm2 suffers the overflow under certain input conditions, and then you may need to pre-condition the input vectors prior to passing on to norm2.

Jim

[Anthony_Richards]

Sorry for my mistake in omitting the square root in the final ACOS statement.
Using the square of the modulus ( and not taking the  square root) for the tests and
putting the square root in the ACOS argument avoids taking an extra square root earlier, thus saving several machine cycles.

[Simon_Geard]

My experience of this in legacy code that gets called millions of times is that atan2 is always better than acos for several reasons:

1. You might think that if you have two normalized vectors _s_ & _t_ then acos(_s_._t_) will always give you the angle between them but you'd be wrong. In practice rounding error can push the dot product outside the valid range so that |s.t| > 1. When this happens the failure is silent - you get an angle but not what you expect. The solution to this is to add guard code, which should be packaged into a function if you're doing it in more than one place, at which point you realize that there ought to be a simpler way and you'd be right - atan2.
2. atan2 is faster than acos (in my tests anyway).
3. I trust that the compiler writers will do a better job than me in handling all the edge cases, accuracy and performance issues which can arise with this calculation and so are best handled by atan2.

[jimdempseyatthecove]

With masked move, X = MIN(X, 1.0), is relatively fast and has no branch instructions. Note, if you already have a test for identity, you might want to use

pragma REAL :: AlmostOne = 1.0 - TINY(1.0)

Tthen use X = MIN(X, AlmostOne)

Of course as to which to use and/or when is subject to your design goal.

Jim Dempsey