Welcome, guest | Sign In | My Account | Store | Cart

Rational approximations (Python recipe) by Raymond Hettinger
ActiveState Code (http://code.activestate.com/recipes/542190/)

Find integer ratio approximations to floating point numbers.

Python, 35 lines
Download
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def approx1(c, maxd):
    'Slow way using sequential search'
    d = min(xrange(2, maxd), key=lambda d: abs(round(c*d)/d/c-1))
    return int(round(c*d)), d

import fractions              # Available in Py2.6 and Py3.0

def approx2(c, maxd):
    'Fast way using continued fractions'
    return fractions.Fraction.from_float(c).limit_denominator(maxd)



if __name__ == '__main__':
    import math
    for x in (math.pi, math.e):
        for d in (10, 100, 1000, 10000, 100000, 1000000):
            print approx1(x, d), approx2(x, d)
        print

##### Sample output
        
(22, 7) 22/7
(311, 99) 311/99
(355, 113) 355/113
(355, 113) 355/113
(312689, 99532) 312689/99532
(3126535, 995207) 3126535/995207
        
(19, 7) 19/7
(193, 71) 193/71
(1457, 536) 1457/536
(25946, 9545) 25946/9545
(271801, 99990) 271801/99990
(1084483, 398959) 1084483/398959

355/113 is the best approximation to Pi with a denominator below 10000. 1457/536 is the best approximation to E with a denominator below 1000.

Tags: algorithms

2 comments

Eduardo Moguillansky 16 years, 3 months ago  # | flag

pi should probably be c

Alexey Marinichev 16 years, 1 month ago  # | flag

Continued fractions work better for this, wikipedia has an article on those. Using continued fractions you can always restore the rational number given its floating point approximation, if its denominator is not too big. Just truncate the fraction if its term gets unreasonably large.

Below is the implementation that is almost always correct, except for those special cases when you must not accept exact half of the last term.

# Generate continued fraction for x.
def chain(x):
    while True:
        a = int(x)
        x = 1/(x-a)
        yield a

# Get best rational approximations of a continued fraction.
def eval_chain(f):
    n1, n2 = 1, 0
    d1, d2 = 0, 1
    for a in f:
        for b in xrange((a+1)/2, a+1):
            yield (b*n1+n2, b*d1+d2)
        n1, n2 = a*n1+n2, n1
        d1, d2 = a*d1+d2, d1

def approx(c, maxd):
    n2, d2 = int(c), 1
    for n1, d1 in eval_chain(chain(c)):
        if d1 > maxd:
            return n2, d2
        n2, d2 = n1, d1
Created by Raymond Hettinger on Wed, 9 Jan 2008 (PSF)
◄ Python recipes (4591) ►
◄ Raymond Hettinger's recipes (97) ►

Required Modules

  • (none specified)

Other Information and Tasks