Best rational approximations of a float

By Gribouillis on Sep 18th, 2009 1:04 pm

This snippet generates the best rational approximations of a floating point number, using the method of continued fractions (tested with python 2.6 and 3.1).

#!/usr/bin/env python

"""A generator of the best rational approximations of a float
   using continued fractions.
   (see http://en.wikipedia.org/wiki/Continued_fraction)

   possible improvements left to the interested reader:
      - extend the generator to support decimal.Decimal numbers.
      - extend the generator to support gmpy.mpf numbers.
"""

def best_rationals(afloat):
  """generate triples (num, den, error) where num/den
  is a best rational approximation of the float afloat and
  error is the difference (afloat - num/den)"""
  afloat, lastnum, num = ((-afloat, -1, int(-afloat)) 
                         if afloat < 0 else (afloat, 1, int(afloat)))
  lastden, den = 0, 1
  rest, quot = afloat, int(afloat)
  while True:
    yield num, den, afloat - float(num)/den
    rest = 1.0/(rest - quot)
    quot = int(rest)
    lastnum, num, lastden, den = (num, quot * num + lastnum,
                                  den, quot * den + lastden)


def testit():
    """Example use of best_rationals()"""
    from math import pi
    value = (pi ** 2) / 6
    gen = best_rationals(value)
    for i in range(20):
        num, den, err = next(gen)
        print ("{0} {1} {2:.2e}" .format(num, den, abs(err)))

if __name__ == "__main__":
    testit()

2 Years Ago

Note that starting with python 2.6, you can use the fractions module, which gives similar results

from math import pi
value = (pi ** 2) / 6
from fractions import Fraction
print(Fraction.from_float(value).limit_denominator(10000000))
"""my output --->
13070401/7945851
"""

Gribouillis

Posting Maven

Moderator

2,786 posts since Jul 2008

Reputation Points: 1,044

Solved Threads: 691

2 Years Ago

If you want approximations for very big numbers, you can use the clnum big numbers library . Here is an exemple showing the square root of 2 with 1000 digits and 2 rational approximations computed with clnum

>>> from clnum import *
>>> s = mpf("2.0", 1000)
>>> v = sqrt(s)
>>> v
mpf('1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372352885092648612494977154218334204285686060146824720771435854874155657069677653720226485447015858801620758474922657226002085584466521458398893944370926591800311388246468157082630100594858704003186480342194897278290641045072636881313739855256117322040245091227700226941127573627280495738108967504018369868368450725799364729060762996941380475654823728997180326802474420629269124859052181004459842150591120249441341728531478105803603371077309182869314710171111683916581726889419758716582152128229518488472089694633862891562883',1019)
>>> ratapx(v, 10000000000000000000000000)
mpq(9165691521498228451812099,6481122629115441680520770)
>>> ratapx(v, 1000000000000000000000000000000000000000000000000)
mpq(848606505733840948241852054437801506505048644099,600055414763419758512382581064986635475955641470)

Gribouillis

Posting Maven

Moderator

2,786 posts since Jul 2008

Reputation Points: 1,044

Solved Threads: 691

2 Years Ago

On the other hand useful fact is also to remember that

>>> print(Fraction.from_float(pi).limit_denominator(10000)) 355/113 >>> 355/113.0 3.1415929203539825 >>> 355/113.0-pi 2.667641894049666e-07 >>>

pyTony

pyMod

Moderator

5,359 posts since Apr 2010

Reputation Points: 782

Solved Threads: 852

1 Year Ago

A nice library to consider is the bigfloat module, which is a python wrapper around the well known mpfr C library based on GNU gmp. Install it from pypi ( easy_install bigfloat ). Then

>>> import bigfloat as bf
>>> with bf.Context(precision=1000):
...  x = bf.BigFloat("2")
...  root = bf.sqrt(x)
...  print root.as_integer_ratio()
... 
(473544376400726394228831039451913480123688799751094630616578363333661533455686216226143010758906440540743599383075282477937325116004184859362870439857750473316274124854227871716121348224479943162567953976022762080355598781683844459583873884188875089875732987543720916348829425487472888348035796043929L, 334846439745708537796382827831250565800439003657979252326171996365734703476542538279124493379904955664873335286735358382870982901778848138624518049209330462622955242963257218294408581408199098183686068192282702343236935664606211486223923248314908216080349889927704442739388432239144512088662677127168L)

Gribouillis

Posting Maven

Moderator

2,786 posts since Jul 2008

Reputation Points: 1,044

Solved Threads: 691

This article has been dead for over three months

You