CODEWITHSUNDEEP
cmath
Michael Sazonov
Michael Sazonov

Reputation: 1533

-1.#IND00 as result of positive root

Here's the code:

double angle( double *vec1 , double *vec2 ){
    return acos( ( *vec1 * *vec2 + *(vec1+1) * *(vec2+1) + *(vec1+2) * *(vec2+2) ) / ( sqrt( pow( *vec1 , 2 ) + pow( *(vec1+1) , 2 ) + pow( *(vec1+2) , 2 ) ) * sqrt( pow( *vec2 , 2 ) + pow( *(vec2+1) , 2 ) + pow( *(vec2+2) , 2 ) ) ) );
}

Vectors entered: <1,1,1> , <2,2,2>

It returns angle: -1.#IND00

Could you please tell me, what's wrong?

Here's the code "readable": 

arcCos[ (x1*x2 + y1*y2 + z1*z2) / sqrt(x1^2 + y1^2 + z1^2) * sqrt(x2^2 + y2^2 + z2^2) ]

Upvotes: 1

Views: 244

Answers (2)

Martin Kristiansen
Martin Kristiansen

Reputation: 10222

ok... here it goes .. the problem is due to math :-)

first I refactored your code :

double angle( double *vec1 , double *vec2 ){

  // sqlen = 6
  double sqlen =  vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]; 

  // sqpow1 = 3
  double sqpow1 = pow( vec1[0] , 2 ) + pow( vec1[1] , 2 ) + pow( vec1[2] , 2 );

  // sqpow2 = 12
  double sqpow2 = pow( vec2[0] , 2 ) + pow( vec2[1] , 2 ) + pow( vec2[2] , 2 );

  // 6/(sqrt(3)*sqrt(12)) = 1
  return acos( sqlen / ( sqrt(sqpow1) * sqrt(sqpow2) ) );
}

This should return '0', since acos(1) is clearly defined as 0. ;-)

Upvotes: 1

Chris
Chris

Reputation: 27619

Your calculation should result in working out acos(1) with the numbers you have given which should come out as 0. However, 1 is at the very top of the allowed range for acos(x) to give a real result, anything higher would be complex (or in most cases throw some kind of error).

My guess would be that you are getting a floating point number slightly greater than one which is causing the problem. The reason is everybody's friend the floating point rounding errors.

When using your method of calculations you will be ending up with the dividend being the product of two square roots. In this case the square roots will of be sqrt(3) and sqrt(12). When multiplied these will give a whole number - 6. However, if both these numbers are rounded down slightly to make them fit in a floating point then when they are multiplied they might be slightly smaller than 6. You will then end up with 6 (on the top) dividided by something slightly smaller than 6 which would give you something slightly bigger than one.

One way to try to mitigate this risk is to change your calculation to:

sqrt((x1^2 + y1^2 + z1^2) * (x2^2 + y2^2 + z2^2))

This has the effect that the sqrt is only being done once on the product which is more likely to produce correct results.

Upvotes: 2

Related Questions