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

[Chris Metzler]

I'm getting some extremely odd results using the trig functions
in the math module.  I don't know whether there's a bug here, or
(more likely) I'm just missing something.  They kinda look like
problems I'd see with truncation/roundoff errors, but they're
occurring for larger values of the floating point exponent than
I'd expect this to happen.

The script below illustrates the kind of thing I'm seeing.  Do
other people see the same thing?  Am I missing something here?

Thanks,

-c



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)

# note that with this formula, as long as x and y are between 0 and pi,
# it's mathematically impossible for gamma1 to be smaller than | x - y |,
# since (1-cos(z)) will range from 0 to 2, and the sines will be positive
# as long as x and y are between 0 and pi.  So I'm always subtracting
# a positive number off of cos(x-y).  So that means the expression
# gives a number smaller than cos(x-y); the arccos should then give an
# angle larger than | x - y |.
# So gamma1 should be greater than | x - y |.  And it is.

print "This should be positive:  ", gamma1 - math.fabs(x-y)

# I get 1.20693498327e-12

# But now let's rearrange the equation.
cos_gamma2 = math.cos(x)*math.cos(y) + math.sin(x)*math.sin(y)*math.cos(z)
gamma2 = math.acos(cos_gamma2)

# This is the SAME EQUATION as above -- simply substitute the
# trigonometric identity cos(x-y) = cos(x)*cos(y) + sin(x)*sin(y)
# into the first equation.
# So, since it's still the same equation, we should get the same
# result -- and in particular, gamma2 should be greater than | x - y |.
# It isn't.

print "This should be positive:  ", gamma2 - math.fabs(x-y)

# I get -3.11013937044e-13

# Note that we don't need anything as weird as my data above
# to generate weirdness.

x = 0.999999
y = 0.999997
z = 1e-08

cos_gamma3 = math.cos(x - y) - math.sin(x)*math.sin(y)*(1.0 - math.cos(z))
gamma3 = math.acos(cos_gamma3)

# Again, just from the functional form, gamma3 should be larger
# than | x - y |.  However . . .

print "This should be positive:  ", gamma3 - math.fabs(x-y)

# it isn't.  I get -2.21217876213e-11.


-- 
Remove the "snip-me" to email.