Created
March 16, 2019 18:59
-
-
Save mgritter/7431977db4cb7a58b9ba44422d3ebe59 to your computer and use it in GitHub Desktop.
Cover's probability test for irrationality
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment