Best rational approximations of a float

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

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

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

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

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

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.