Entry

module for continued fractions, useful for finding optimal fractional approximation to real numbers

Jul 11th, 2000 04:19
unknown unknown, Huaiyu Zhu

#!/usr/bin/env python
# (C) 2000 Huaiyu Zhu <hzhu@users.sourceforge.net>.      Licence: GPL
# $Id: ContFrac.py,v 1.3 2021/07/10 05:09:11 hzhu Exp $
"""
ContFrac.py - Continued fraction approximation utilities
        Convert between float, continued fraction and fractions.
        Useful for finding optimal rational approximation at any 
precision.
Usage examples: python ContFrac.py
Representation:
        Continued fraction: list of integers.
        Ordinary fraction: pair of long integers. 
Optional variables in the functions:
        nterm: max number of terms in continued fraction.
        prec:  precition (inverse relative error) - maybe too 
conservative?
"""
def cf2real(cf):
        "Convert to float from continued fraction"
        x = 0.
        for i in range(1,len(cf)+1):            x = 1./(cf[-i] + x)
        return  1./x
def real2cf(x, prec=1e16, nterms=40):
        "Convert to continued fraction from float"
        cf = [];        p = 1.; a = 1
        for i in range(nterms):
                a, a0 = int(x), long(a)
                cf.append(a)
                p = p * (abs(a*a0) + 1)
                if p > prec or abs(x-a)<1e-12: break
                x = 1./(x-a)
        "Merge last term if it is one, so the expression is unique"
        n = len(cf)-1
        if cf[n] == 1: cf[n-1]=cf[n-1]+1; del cf[n]
        return cf
def cf2frc(cf, prec=1e16):
        "Convert to fraction from continued fraction"
        "Estimates term n that achieves precision - need better 
algorithm"
        p = 1.; n = 0
        for n in range(1, len(cf)):
                p = p * (cf[n-1]*(cf[n]+.5)+1)
                if p >= prec: break
        "Conversion using first n terms"
        x, y = long(cf[n]), 1L
        for i in range(n):
                x, y = x*cf[n-1-i] + y, x
        return  x, y
def frc2cf(frc):
        "Convert to continued fraction from fraction"
        cf = []
        x, y = frc
        while 1:
                if y == 0: break
                a, b  = x/y, x%y
                cf.append(a)
                x, y = y, b
        return cf
#------------------------------------------------------------------
def test(x, prec, nterms):
        print "Convert %.16f to continued fraction and back to float" % 
x
        cf = real2cf(x, nterms=nterms)
        print cf,  cf2real(cf)
        print "Convert to fraction and back to float, limit precision: 
%g" %prec
        print "num. of terms                     fraction       float   
        error"
        for n in range(1, nterms+1):
                f = cf2frc(cf[:n], prec)
                r = f[0]/float(f[1])
                print "%2d      %33s  % .16f  %- 8.4g" % (n, f, r, 
r/x-1)
if __name__ == "__main__":
        from math import pi, exp, sqrt
        print "------------ Effect of too much precision and terms 
------------"
        test(-0.6667, prec=1e19, nterms=6)
        print "------------ Effect of limited number of terms 
----------------"
        test(pi, prec=1e16, nterms=6)
        print "------------ Effect of limited of precision 
--------------------"
        test(exp(1), prec=1e4, nterms=9)
        # Interesting pattern in continued fraction of e.
        print "------------ Conversion between fraction and continued 
fraction"
        cf = real2cf(1/sqrt(2), nterms=17); print cf
        f = cf2frc(cf); print f
        cf1 = frc2cf(f); print cf1
        if cf == cf1: print "Idential"
        else: print "Something wrong"