Arccos, out of range problem

Hello (sorry for my poor english)

I need to calculate the angle theta between vectors v1 and v2 using the scalar product relation :

theta=acos( v1*v2/( norm(v1)*norm(v2) ) )

So I implemented the relation using a three double x,y,z structure for the vectors (physical sens) and i got the following output.

Nv1: 1, Nv2: 1, prod: -1    
theta: nan


Where Nv1,2 means norm(v1,2) and prod means v1*v2 the scalar product.

I guess the error comes from the precision of the double. I mean during the execution -1/1*1 = -1 - some remaining value in the last digit.

Do you know an easy way the fix it? I do this operation in a large set of uniformly randomly rotated units vectors in 3D.

Thank you for your help.
Can you show us your code?
Arccos out of range?

Well, the acceptable range of inputs is: -1 <= value <= +1
could you insert some test prior to calling arccos, to validate the input.
e.g. if (value < -1) value = -1; and the same for +1.
Last edited on
Thank you for your reactions,

I'll use "if" statement to ensure that the value remains between -1 and 1.

I was just wondering if there was a better way to avoid numerical inaccuracy.



It would be good to see some of your code, so we can see if there are any problems with your methodology.
My code looks like :
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <math.h>

struct vec 
{
	double x;
	double y;
	double z;
};

double get_theta(vec v, vec w){
     double Nv=sqrt( pow(v.x,2) + pow(v.y,2) + pow(v.z,2));
     double Nw=sqrt( pow(w.x,2) + pow(w.y,2) + pow(w.z,2));
     double Scal_prod=v.x*w.x+v.y*w.y+v.z*w.z;
     cout<<"Nv1: "<<Nv<<" , Nv2: "<<Nw<<", prod: "<<prod<<endl;
     returun acos( Scal_prod/( Nv*Nw ) );
} 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// In main e.g.

vec v,w;

v.x=((double)rand()/(double)RAND_MAX);
v.y=((double)rand()/(double)RAND_MAX);
v.z=sqrt(1-pow(v.x,2)-pow(v.y,2));

w.x=((double)rand()/(double)RAND_MAX);
w.y=((double)rand()/(double)RAND_MAX);
w.z=sqrt(1-pow(w.x,2)-pow(w.y,2));

theta=get_theta(v,w);
cout<<"theta: "<<theta;


In this example I don't need to compute the norms (=1), but I use this function in more general case.

For some value of v and w, I get the output I showed in my first post.

(Mathematically, it's easy to prove that for all v,w in R3, v*w/(norm(v)norm(w)) is in [-1,1] for every norm built from a scalar product

by the way, can LaTex be used to write equations in post?)

Just some terminology that seems different from when I learnt vectors at Uni.

In my mind, This:

double Nv=sqrt( pow(v.x,2) + pow(v.y,2) + pow(v.z,2));

calculates the magnitude or length of vector v.

and this:
double Scal_prod=v.x*w.x+v.y*w.y+v.z*w.z;

is a dot product, not a scalar product. Scalar product is when you multiply each member of the vector by one number - the scalar, it has the effect of "scaling up" the vector.

Never mind, I have just read the wiki page. However I still think my descriptions of scalar, dot, magnitude, unit vector & normal are easier to understand.

I have heard people use the term "normalise", but from context I took it to mean unit vector or otherwise a magnitude. I find this confusing because it is too near to the meaning of a "normal vector", which means "at right angles to" as in a normal to a plane.

Does your naming of the variable Nv refer to the normal of vector v? It is the magnitude of v at the moment.

after all that, to your problem:

The variable prod is not initialised to anything in the function get_theta - did you mean Scal_prod? Probably just a typo.

I can't see anything really wrong with your code, so you should try firstly to use known values rather than random ones.

Next, use the debugger with a watchlist of variables to see where it goes wrong - this is by far the easiest way.

Failing that, put loads of cout statements to check on the variable values.

Interested to see how you go. Cheers

Edit: when I wrote some code to deal with vector I had functions that calculated scalar, dot, magnitude, unit vector & normal, because these functions are used all the time with vector operations.

Edit2:

This is from the reference page:

If the argument is out of this interval, a domain error occurs, setting the global variable errno to the value EDOM.


Can you make use of this to see if the acos function worked?
Last edited on
In my experience you can get these domain errors when the argument is a tiny smidge outside the [-1,1] interval.
Usually when the vectors involved should be unit vectors but have mutated a bit from some long calculation.
The only way is to force these back into [-1,1], I think.

why not Nv=sqrt(v.x*v.x+v.y*v.y+v.z*v.z) instead of using pow ?
Topic archived. No new replies allowed.