Some notes about float approximations in mxNumber

[Kay Schluehr]

Hi Mark,

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