Floating-point glitches with the math module. Bug? Or am I missing something?

[Alexandre Fayolle]

Le 20-09-2004, Chris <cmetzler at speakeasy.snip-me.net> nous disait:
> import math
>
> x = 0.90455926850657064620
> y = 0.90463240818925083619
> z = -2.00200807043415807129e-08
>
>
> cos_gamma1 = math.cos(x - y) - math.sin(x)*math.sin(y)*(1.0 - math.cos(z))
> gamma1 = math.acos(cos_gamma1)

Here's my attempt. 

Firstly, I use the fact that scipy provides a cosm1 function which 
returns cos(x)-1, and which is more precise that for small values of x. 

>>> from scipy import special as S
>>> cos_gamma_minus_1 = S.cosm1(x-y) - S.sin(x)*S.sin-y)*S.cosm1(z)
>>> print cos_gamma_minus_1
-2.67470646622e-09

Then I'd like to use an acosm1 function in order to keep the precision.
Since there is no such function available, I'll just solve
cosm1(x)-cos_gamma_minus_1=0 for which I know there is a solution
between 0 and 1. Any solver will do the trick. To keep things simple,
I'll use scipy.optimize.zeros.bisect (not the fastest around, but does a
good job)

>>> from scipy.optimize import zeros
>>> gamma = zeros.bisect(lambda x: S.cosm1(x)-cos_gamma1_minus1, 
...                      0., 1., 1e-20)
>>> print gamma
7.31396809867e-05

-- 
Alexandre Fayolle                              LOGILAB, Paris (France).
http://www.logilab.com   http://www.logilab.fr  http://www.logilab.org