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

[Tim Peters]

[Jeff Epler]
> ...
> acos(cos(d)) for small d gives fairly large errors on my system:
> >>> d = .999999 - .9999977
> >>> acos(cos(d)), d
> (1.2999928800133351e-06,
>  1.2999999999818712e-06)
> >>> #   ^ First difference here
> a half dozen correct digits isn't great, but it's all my platform's trig
> functions seem to give me.

This is an instance of a general phenomenon.  Compute
cos(acos(cos(d))), and compare it to cos(d).  On my box:

>>> cos(acos(cos(d))) == cos(d)
True
>>>

But I also have:

>>> d
1.2999999999818712e-006
>>> acos(cos(d))
1.2999928800130605e-006
>>>

But as shown by the first line, cos() applied to either of those
computes exactly the same result -- to machine precision, d and
acos(cos(d)) have equal claim to being "the correct" inverse of
cos(d).

This isn't the fault of sloppy trig functions, it's an unavoidable
consequence of using finite precision arithmetic to evaluate a
function *near* a local minimum or local maximum.  cos() has a local
maximum at cos(0), and d is near 0.  For any reasonable function f
evaluated at point a such that f(a) is a local min or max, the first
derivative of f at a is 0.  Therefore f(a+h) ~= f(a) + h**2*f''(a)/2
(the first-order term of the Taylor expansion around a vanishes). 
That in turn roughly means that if you change any of the bits in the
least significant half of a, it makes no visible difference to the
computed value of f(a).

That's why any numerical analysis text will tell you that you can't
expect to compute the value at which a function achieves a local min
or max to better than about half the bits of precision in your
floating-point format.  Computing inverses near such peaks/valleys
suffers the same problem, for the same reason:  while the mathematical
cos() is 1-to-1 near 0, the machine-precision cos() is many-to-1 near
0, because cos() is very flat near 0 (== the first derivative is 0 at
cos(0)).  And since machine-precision cos() is many-to-1 near 0,
machine-precision acos() near cos(d) for a small |d| can pick many
results e for which cos(e) is exactly equal to cos(d) to machine
precision.

Note that your box gave a different acos(cos(d)) than my box did.  Nevertheless,

>>> cos(d) == cos(acos(cos(d))) == cos(1.2999928800133351e-06)
True
>>>

on my box, where the far right is the answer your box gave for
acos(cos(d)).  There are simply a ton of floats that have the same
machine-precision cosine.