Arccos, out of range problem

[Tim1294]

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.

[TheIdeasMan]

Can you show us your code?

[Chervil]

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.

[Tim1294]

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.

[TheIdeasMan]

It would be good to see some of your code, so we can see if there are any problems with your methodology.

[Tim1294]

My code looks like :

#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 ) );
} 

// 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 R³, 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?)

[TheIdeasMan]

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?

[mik2718]

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 ?