Some notes about float approximations in mxNumber

[M.-A. Lemburg]

Kay Schluehr wrote:
> Hi Marc,
> 
> I was a bit surprised to find the very slow Farey approximation by
> means of the <FareyRational> class in the mxNumber package. If the goal
> was to reconstruct a rational from a float it is not a good choice and
> should be replaced by a continued fractions approximation. 

The idea was to be able to create a Rational() from a float
using a given upper bound on the denominator.

The standard Rational() constructor already provides a way
to construct a rational out of a Python or mxNumber float,
but always uses maximum precision - which may not always be
what the user wants (e.g. to work around rounding errors).

Thanks for posting your faster version. I think it would be
a good candidate for a Python Cookbook Entry and I'll
see whether I can add something like this to one of the next
mxNumber releases.

> Some time
> ago I implemented it by myself so it can be published here:
> 
> 
> def cfrac(z,bound=10**-5,n=10):
>     '''
>     creates a continued fraction from some number z.
>     @ bound - terminate cf if the rest is lower than the provided bound
>     @ n     - terminate cf after n steps
>     '''
>     l = []
>     while 1:
>         a = int(z)
>         l.append(a)
>         y = z-a
>         if y<=bound or len(l)==n:
>             return l
>         z = 1./y
> 
> 
> def fold(cf):
>     '''
>     create u/v = cfrac(a0,a1,...,an) using following rules:
> 
>     1. / u1 u0 \      / a1  1 \    / a0  1 \
>        |       |  =   |       | *  |       |
>        \ v1 v0 /      \ 1   0 /    \ 1   0 /
> 
>     2. The recursion rules
>        v(n+1) = v(n)*a(n+1)+v(n-1)
>        u(n+1) = u(n)*a(n+1)+u(n-1)
>     '''
>     if len(cf)<2:
>         return Rational(0,1)
>     un   = cf[0]*cf[1]+1
>     vn   = cf[1]
>     un_1 = cf[0]
>     vn_1 = 1
>     for a in cf[2:]:
>         b  = un
>         un = un*a+un_1
>         un_1 = b
>         b  = vn
>         vn = vn*a+vn_1
>         vn_1 = b
>     return Rational(un,vn)
> 
> 
> 
>>>>fold(cfrac(1./3))
> 
> 1/3
> 
>>>>fold(cfract(1525/42.))
> 
> 1525/42
> 
> import math
> 
> 
>>>>fold(cfract(math.sqrt(2)))
> 
> 3363/2378
> 
> 
> It is possible to provide this functionality in an even more efficient
> manner because it is usefull to bound only the approximation error of
> the rational, not the size of the continued fraction.
> 
> def float2ratio(z, bound=10**-5):
>     '''
>     convert a float into a Rational.
>     The approximation is bound by the error-limit 'bound'.
>     '''
>     a = int(z)
>     y = z-a
>     z = 1./y
>     b = int(z)
>     un   = a*b+1
>     vn   = b
>     un_1 = a
>     vn_1 = 1
>     a    = b
>     while 1:
>         y = z-a
>         if y<bound:
>             return Rational(un,vn),k
>         z = 1./y
>         a = int(z)
>         x  = un
>         un = un*a+un_1
>         un_1 = x
>         x  = vn
>         vn = vn*a+vn_1
>         vn_1 = x
>         xn = float(un)/vn
>         yn = float(un_1)/vn_1
>         if abs(xn-yn)<=bound:
>             return Rational(un,vn)
> 
> 
>>>>math.sqrt(2)
> 
> 1.4142135623730951
> 
> 
>>>>float2ratio(math.sqrt(2))
> 
> 1393/985
> 
>>>>1393./985
> 
> 1.4142131979695431
>         ^
> 
> 
>>>>float2ratio(math.sqrt(2),10**-10)
> 
> 275807/195025
> 
>>>>275807./195025
> 
> 1.4142135623637995
>             ^
> 
> 
>>>>math.pi
> 
> 3.1415926535897931
> 
> 
>>>>float2ratio(math.pi,10**-14)
> 
> 245850922/78256779
> 
>>>>245850922./78256779
> 
> 3.1415926535897931
> 
> Note that the algorithm needed only 13 iteration steps to approximate
> pi 
> in this accuracy.
> 
> Regards,
> Kay
> 

-- 
Marc-Andre Lemburg
eGenix.com

Professional Python Services directly from the Source  (#1, Apr 04 2005)
>>> Python/Zope Consulting and Support ...        http://www.egenix.com/
>>> mxODBC.Zope.Database.Adapter ...             http://zope.egenix.com/
>>> mxODBC, mxDateTime, mxTextTools ...        http://python.egenix.com/
________________________________________________________________________

::: Try mxODBC.Zope.DA for Windows,Linux,Solaris,FreeBSD for free ! ::::