Cover's probability test for irrationality

cover.py

import mpmath
from fractions import Fraction

def earliestDecimalMatch( f, nDigits ):
    """Return the earliest fraction in the enumeration
       0/1, 1/1, 1/2, 1/3, 2/3, 1/4, 2/4, 3/4, 1/5, ...
       that matches the value of f to nDigits of accuracy."""

    # 3 digits: 0.333 * 1000 = 333
    # Clueless search: try every possible denominator from smallest to largest
    maxDenom = 10 ** nDigits
    for d in xrange(1, maxDenom + 1):
        # Consider 0.951 with n=1.  the nearest integer is 10, but 10/10
        # doesn't match.  9/10 does.
        # We might have to go for the ceiling sometimes, though???
        n = mpmath.floor( mpmath.fmul( f, d ) )

        candidate = mpmath.floor( mpmath.fmul( mpmath.fdiv( n, d ),
                                               maxDenom ) )
        actual = mpmath.floor( mpmath.fmul( f, maxDenom ) )
        #print d, n, candidate, actual
        
        if candidate == actual:
            return (int(n),d)

        n = n + 1
        candidate = mpmath.floor( mpmath.fmul( mpmath.fdiv( n, d ),
                                               maxDenom ) )
        #print d, n, candidate, actual
        if candidate == actual:
            return (int(n),d)
    return (1,1)

def earliestDecimalMatch_zigzag( f, nDigits ):
    """Return the earliest fraction in the enumeration
       0/1, 1/1, 1/2, 1/3, 2/3, 1/4, 2/4, 3/4, 1/5, ...
       that matches the value of f to nDigits of accuracy."""

    # Start at denominator 0/1.  If number is too small, increase numerator.
    # If number too larger, increase denominator.
    maxDenom = 10 ** nDigits
    d = 1
    n = 0
    actual = mpmath.floor( mpmath.fmul( f, maxDenom ) )
    while d <= maxDenom:
        candidate = mpmath.floor( mpmath.fraction( n * maxDenom, d ) )
        #print d, n, candidate, actual
        if candidate == actual:
            return (n,d)
        elif candidate < actual:
            n += 1
        else:
            d += 1
    return (1,1)

def continuedFraction( r ):
    """Given a Rational number "r", return the continued fraction
    representation as a list of integers."""

    integerPart = r.numerator // r.denominator
    fractionalPart = r - integerPart
    if fractionalPart == 0:
        return [integerPart]
    else:
        return [integerPart] + continuedFraction( 1 / fractionalPart )

def continuedFractionToRational( cf ):
    """Given a finite continued fraction [a_0; a_1, a_2, ..., a_n]
    return it as a rational number."""
    if len( cf ) == 1:
        return Fraction( cf[0], 1 )

    return cf[0] + Fraction( 1, continuedFractionToRational( cf[1:] ) )

def bestRationalInInterval( r1, r2 ):
    """Given two rational numbers, return the best rational, i.e., smallest
    denominator and smallest numerator, between the two."""
    cf1 = continuedFraction( r1 )
    cf2 = continuedFraction( r2 )
    result = []
    # Append "plus infinity": to the shorter continued fraction
    if len( cf1 ) < len( cf2 ):
        cf1 += [ cf2[-1] + 100 ]
    elif len( cf1 ) > len( cf2 ):
        cf2 += [ cf1[-1] + 100 ]
    #print cf1
    #print cf2
        
    for n1, n2 in zip( cf1, cf2 ):
        if n1 == n2:
            result.append( n1 )
        else:
            result.append( min( n1, n2 ) + 1 )
            break
        
    #print result
    between = continuedFractionToRational( result )
    if between.denominator > r1.denominator:
        return r1
    else:
        return between

def earliestDecimalMatch_cf( f, nDigits ):
    """Return the earliest fraction in the enumeration
       0/1, 1/1, 1/2, 1/3, 2/3, 1/4, 2/4, 3/4, 1/5, ...
       that matches the value of f to nDigits of accuracy.""" 
    maxDenom = 10 ** nDigits
    # for example with f = 0.31415926535
    # lower = loor ( f * 1000 ) = 314
    # look for the best rational between 314/1000 and 315/1000
    # This method doesn't pick the endpoints, but sometimes it should.
    lower = mpmath.floor( mpmath.fmul( f, maxDenom ) )
    upper = lower + 1
    frac = bestRationalInInterval( Fraction( int( lower ), maxDenom ),
                                   Fraction( int( upper ), maxDenom ) )
    return (frac.numerator, frac.denominator)
    
def test():
    print "Fractions close to pi/10"
    f = mpmath.mpf( "0.3141592653589" )
    print earliestDecimalMatch( f, 1 )
    print earliestDecimalMatch( f, 2 )
    print earliestDecimalMatch( f, 3 )
    print earliestDecimalMatch( f, 4 )
    print earliestDecimalMatch( f, 5 )
    print earliestDecimalMatch( f, 6 )

    print "Fractions close to 13/1117"
    f = mpmath.fraction( 13, 1117 )
    print earliestDecimalMatch( f, 1 )
    print earliestDecimalMatch( f, 2 )
    print earliestDecimalMatch( f, 3 )
    print earliestDecimalMatch( f, 4 )
    print earliestDecimalMatch( f, 5 )
    print earliestDecimalMatch( f, 6 )    

def indexOfFraction( n, d ):
    """Return the index of the fraction n/d in the enumeration
    0/1, 1/1, 1/2, 1/3, 2/3, 1/4, 2/4, 3/4, 1/5, ...
    """

    if n == 0:
        return 1
    if d == 1:
        return 2
    
    # d=1 is index 2--2
    # d=2 is index 3--3   (+1)
    # d=3 is index 4--5   (+1)
    # d=4 is index 6--8   (+2)
    # d=5 is index 9--12  (+3)
    # ...
    # start(d) = start(d-1) + (d-2)
    # start(2) = 3
    # start(3) = 3 + 1
    # start(4) = 3 + 1 + 2    
    # start(d) = 3 + sum_{i=1}^{i=d-2} i
    #          = 3 + (d-2)(d-1)/2
    # start(4) = 3 + 2*3/2 = 6
    # start(6) = 3 + 4*5/2 = 13

    start = 3 + (d-2)*(d-1) // 2
    return start + (n-1)

def covers_test( f, n, verbose=True ):
    (s_n, s_d) = earliestDecimalMatch_cf( f, n )
    i_n = indexOfFraction( s_n, s_d )
    k_n = 9 ** n
    if verbose:
        # mpmath.nstr rounds instead of truncating.
        truncated = mpmath.fdiv( mpmath.floor( mpmath.fmul( f, 10 ** n ) ),
                                 10 ** n )
        print "f={} s_n={}/{}, i_n={} {}".format(
                mpmath.nstr( truncated, n ),
                s_n, s_d,
                i_n,
                i_n <= k_n )
        
    return i_n <= k_n

def reproduce_paper_results():
    mpmath.mp.dps = 11

    e_over_10 = mpmath.fdiv( mpmath.exp(1), 10 )
    covers_test( e_over_10, 1, True )
    covers_test( e_over_10, 5, True )
    covers_test( e_over_10, 9, True )
    
    pi_over_10 = mpmath.fdiv( mpmath.pi, 10 )
    covers_test( pi_over_10, 1, True )
    covers_test( pi_over_10, 5, True )
    covers_test( pi_over_10, 9, True )

    one_over_7 = mpmath.fraction( 1, 7 )
    covers_test( one_over_7, 1, True )
    covers_test( one_over_7, 5, True )
    covers_test( one_over_7, 9, True )
    
    
def cosine_test( intervals, maxDigits = 9 ):
    mpmath.mp.dps = maxDigits + 2
    pi_over_2 = mpmath.fdiv( mpmath.pi, 2 )
    headerFormat = "cos( {}*pi/{} )"
    headerLen = len( headerFormat.format( intervals-1, 2*intervals ) )
    print (" " * headerLen), " ".join( "{:2}".format( i ) for i in range( 1, maxDigits + 1 ) )

    for i in range( 1, intervals ):
        theta = mpmath.fmul( pi_over_2, mpmath.fraction( i, intervals ) )
        cos_theta = mpmath.cos( theta )
        header = headerFormat.format( i, 2*intervals )
        print ("{:" + str(headerLen) + "}").format( header ),
        for n in range( 1, maxDigits + 1 ):
            if covers_test( cos_theta, n, False ):
                print " R",
            else:
                print " I",
        print cos_theta
        
        
def rationals_test( maxDigits = 9 ):
    mpmath.mp.dps = maxDigits + 2
    rationals = [
        mpmath.fraction( 1, 137 ),
        mpmath.fraction( 1, 2**13 ),
        mpmath.fraction( 4, 7 ),
        mpmath.fraction( 999999999, 1000000000 ),
        mpmath.fraction( 1, 97 ),
        mpmath.fraction( 4, 61 ),
        mpmath.fraction( 54, 59 )
    ]

    print " ".join( "{:2}".format( i ) for i in range( 1, maxDigits + 1 ) )
    for u in rationals:
        for n in range( 1, maxDigits + 1 ):
            if covers_test( u, n, False ):
                print " R",
            else:
                print " I",
        print u