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

[Tim Peters]

[Tim, sez that when you evaluate a function near a local min or max,
 then it's unavoidable that using finite-precision arithmetic makes
 a 1-to-1 mathematical function appear to be many-to-1]

[Terry Reedy]
> To illustrate:
> >>> for i in range(50): print `cos(i*1e-9)` # Windows, 2.2
> 1.0
> 1.0
> 1.0
> 1.0
> 1.0
> 1.0
> 1.0
> 1.0
> 1.0
> 1.0
> 1.0
> 0.99999999999999989
> 0.99999999999999989

...

> The mathematically smooth cos(x) is effectively a stairstep function for
> |x|<5e-8.  Stairsteps are hard to invert ;-).

Yup!  I'd like to point out too that the new decimal module lowers the
bar for doing numeric *investigation* in two important ways:  (1) just
because it's decimal, it's much easier for people to understand what
they're seeing; and, (2) because you can boost precision to anything
you like whenever you like, it's often dead easy to write a function
that delivers excellent accuracy.

I won't say much about the following, but hope it illustrates both
points for those who run it and stare at it.  The decimal cosine
routine is about as naive as you can imagine, but by temporarily
boosting the working precision to twice the caller's precision, then
rounding back at the end, it delivers high-quality results.  Until the
input argument gets "too big".  Then it's a disaster.  So I hope this
also illustrates that decimal isn't a magical cure either.

"""
import decimal, math

def deccos(x):
    origprec = decimal.getcontext().prec
    decimal.getcontext().prec = 2 * origprec
    result = term = decimal.Decimal(1)
    # Compute 1 - x**2/2! + x**4/4! - x**6/6! ...
    i = 1
    sq = -x*x
    while True:
        term *= sq
        term /= i*(i+1)
        i += 2
        newresult = result + term
        if result == newresult:
            break
        result = newresult
    decimal.getcontext().prec = origprec
    return result + 0 # round to orig precision

decimal.getcontext().prec = 20
base = decimal.Decimal("1e-9")

for i in range(50):
    x = base * i
    dcos = deccos(x)
    fcos = math.cos(float(x))
    print "%6s %22s %19.17f %s" % (x,
                                   dcos,
                                   fcos,
                                   fcos == float(dcos))
"""

One thing to notice is that, working with 20 decimal digits, all these
inputs have distinct computed cos() values.  Another is that
math.cos() delivers the same results, after rounding to native binary
fp precision.  Another is that while we can tell that's so because the
4th column always says "True", you'd have a hard time guessing that
from staring at the 2nd and 3rd columns; for example, it's probably
not intuitively obvious that the decimal

    0.999999999999999838
"rounds to" 
    0.99999999999999989

when converted to native binary fp and then back to 17-digit decimal. 
Getting binary<->decimal conversion errors out of the equation makes
working with the decimal module easier for everyone.

decimal-is-too-important-to-waste-on-money<wink>-ly y'rs  - tim