4

I'm very confused with the issue. The code I am running is:

double dotProduct = dot(A, B);
std::cout << dotProduct << std::endl;
theta = acos(dotProduct);
std::cout << theta << std::endl;

The output of which is

-1
ANGLE: -nan(ind)

However, the following works:

double dotProduct = dot(A, B);
std::cout << dotProduct << std::endl;
dotProduct = -1;
theta = acos(dotProduct);

Giving the output:

-1
ANGLE: 3.14159

Also, if I cast dotProduct to a float, acos() outputs the angle correctly:

double dotProduct = dot(A, B);
std::cout << dotProduct << std::endl;
theta = acos((float) dotProduct);

Also resulting in the output

-1
ANGLE: 3.14159

For the dot() function, I am using the Armadillo library. What I don't understand is why acos() should work when I set dotProduct = -1, but will not when it is outputted by the dot() function. What am I doing wrong?

Brendan M
  • 51
  • 1
  • 4
  • 6
    Likely the value returned by `dot_product` is slightly greater than one, and the default format for inserting a `double` into an output stream does not show the value fully. Inserting `std::setprecision(99)` before the `double` should show it. Converting it to `float` rounds to the nearest representable `float`, which is exactly one because the value is so close to one. – Eric Postpischil Sep 02 '18 at 15:01
  • 2
    If @EricPostpischil is right, you can use `dotProduct = std::clamp(dotProduct, -1., 1.);`. – Evg Sep 02 '18 at 15:07
  • 3
    (Note my “slightly greater than one” should be “slightly greater than one in magnitude.” E.g., −1.00001…) – Eric Postpischil Sep 02 '18 at 15:08

1 Answers1

5

I suppose that A and B are normalized vectors, and therefore you expect that dot(A, B) will be between -1 and 1. With floating-point math, this is not necessarily true. |dot(A, B)| can be slightly larger than 1. I bet, if you print it with larger precision, you'll see that your value is slightly smaller than -1.

So, you'll need to clamp dot(A, B) between -1 and 1. Or even better, call acos() only, if the input is between -1 and 1:

double safe_acos(double value) {
    if (value<=-1.0) {
        return <pi>;
    } else if (value>=1.0) {
        return 0;
    } else {
        return acos(value);
    }
}

This could be slightly faster, as you avoid calling acos() for out-of-bound values.

Note: If you do 3D math, there is a high chance that you can avoid calling acos() altogether. Using trigonometric identities, usually we can replace it with something faster, and more accurate solution.

geza
  • 28,403
  • 6
  • 61
  • 135