[Numpy-discussion] Arcos returns nan

[David Cournapeau]

David Cournapeau wrote:
> João Quinta da Fonseca wrote:
>   
>> I have a function that returns the dot product of two unit vectors.  
>> When I try to use arcos on the values returned I sometimes get the  
>> warning:
>> "Warning: invalid value encountered in arccos", and the angle  
>> returned is nan. I found out that this happens for essentially co- 
>> linear vectors, for which the dot product function returns 1.0.  This  
>> looks like 1.0 no matter how I print it but if I do: >>N.dot(a,b)>1,  
>> I get: >>True.
>> Now I guess this arises because of the way computers store floats but  
>> shouldn't numpy take care of this somehow? Is it a bug? I don't seem  
>> to have this problem with Matlab or Fortran.
>>     

Note that the below program in C has the same behaviour than numpy if I 
understand your situation right:

#include <math.h>
#include <float.h>
#include <stdio.h>

int main()
{      
        fprintf(stderr, "%f\n", acos(1.));
        fprintf(stderr, "%f\n", acos(1. + DBL_EPSILON));

        return 0;
}    

(EPSILON being by definition the smallest value such as 1. + EPSILON > 
1.). So the problem in your case is likely to be caused by precision 
problems when you normalized the scalar product. An easy solution would 
be to clip all the values > 1. to 1 (same for values < -1); another 
solution may be in a better way to do the normalization.

But I don't see how matlab or Fortran would be any different.

cheers,

David