[Python-Dev] Picking on platform fmod

[Tim Peters]

Here's your chance to prove your favorite platform isn't a worthless pile of
monkey crap <wink>.  Please run the attached.  If it prints anything other
than

0 failures in 10000 tries

it will probably print a lot.  In that case I'd like to know which flavor of
C+libc+libm you're using, and the OS; a few of the failures it prints may be
helpful too.  If it only prints one or two failures, it's probably a bug in
*my* code, so I especially want to know about that.  If it dies with an
assertion error, that would be mondo interesting to know, and then the
flavor of HW chip may also be relevant.

I already know MSVC 6 under Win98SE passes ("0 failures") on two distinct
boxes <wink>, so no need for more about that one.

What you're looking for:  Given finite doubles x>0 and y>0, there's a unique
integer N and unique real number R such that

    x = N*y + R

exactly and 0 <= R < y.  N may not be *representable* as an integer (or
long, or long long, or double) on the machine, but it's A Theorem that the
infinitely-precise value of R is exactly representable as a machine double.

The C stds have never (IMO) been clear about whether fmod(x, y) must return
this exact R, although the C99 *rationale* is clear that this is the intent
(why committees don't fold these subtleties into the bodies of their stds is
beyond me).

The program below generates nasty test cases at random, computes the exact R
in a clever (read "on the edge of not working but pretty fast") way using
Python, and compares it to the platform fmod result.  It takes about 6
seconds to run on a 866MHz box using current CVS Python, so it shouldn't be
a burden to try.

If you get no failure, of course I'd like to hear about that too.

it's-not-like-you-had-anything-fun-to-do-this-weekend-ly y'rs  - tim


from math import frexp as _frexp, ldexp as _ldexp

# ffmod is a Pythonic variant of fmod, returning a remainder with the
# same sign as y.  Excepting infs and NaNs, the result is exact.

def ffmod(x, y):
    if y == 0:
        raise ZeroDivisionError("ffmod: divide by 0")
    remainder = abs(x)
    yabs = abs(y)
    if remainder >= yabs:
        dexp = _frexp(remainder)[1] - _frexp(yabs)[1]
        assert dexp >= 0
        yshifted = _ldexp(yabs, dexp)  # exact
        for i in xrange(dexp + 1):
            # compute one bit of the quotient (but not materialized;
            # we only care about the remainder at the end)
            if remainder >= yshifted:
                assert remainder < yshifted * 2.0
                remainder -= yshifted  # exact
            yshifted *= 0.5            # exact
        assert yshifted * 2.0 == yabs
    assert remainder < yabs
    if y < 0 and remainder > 0:
        remainder -= yabs              # exact
        assert remainder < 0
    return remainder

# ffmod and C99 fmod should agree whenever x>0 and y>0.  Try one, and
# return 1 iff they don't agree.

def _tryone(x, y, dump=0):
    n = math.fmod(x, y)
    e = ffmod(x, y)
    if dump:
        print "fmod" + `x, y`
        print "native:", `n`
        print " exact:", `e`
    return n != e

# Test n random inputs, in the sense of random mantissas and random
# exponents in range(-300, 301).  The hardest cases have x much larger
# than y, and this will generate lots of those.

def _test(n, showresults=0):
    from random import random, randrange
    nfail = 0
    for i in xrange(n):
        x = _ldexp(random(), randrange(-300, 301))
        y = _ldexp(random(), randrange(-300, 301))
        if x < y:
            x, y = y, x
        if _tryone(x, y, showresults):
            nfail += 1
            _tryone(x, y, 1)
    print nfail, "failures in", n, "tries"

if __name__ == "__main__":
    import math
    _test(10000, 0)