Best rational approximations of a float

Gribouillis

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).

748 Views
About the Author

Mind explorer.

#!/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()
Gribouillis 1,391 Programming Explorer Team Colleague

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 1,391 Programming Explorer Team Colleague

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)
TrustyTony 888 pyMod Team Colleague Featured Poster

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
>>>

Gribouillis 1,391 Programming Explorer Team Colleague

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)
Be a part of the DaniWeb community

We're a friendly, industry-focused community of 1.19 million developers, IT pros, digital marketers, and technology enthusiasts learning and sharing knowledge.