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 ! :::: -- http://mail.python.org/mailman/listinfo/python-list
