At 10:40 AM 1/6/2007, Dick Moores wrote: >At 11:21 PM 1/4/2007, Terry Carroll wrote: > >On Wed, 3 Jan 2007, Dick Moores wrote: > > > > > Be that as it may, farey() is an amazing program. > > > >Not to beat this subject to death, but the comment at the bottom of > >http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/52317 about > >continued fractions piqued my interest. I'm no mathematician, but I > >encountered continued fractions a long time ago and was fascinated by > >them. So I read the URL pointed to, > >http://mathworld.wolfram.com/ContinuedFraction.html , and came up with the > >following: > > > >##################################################### > > > >def cf(x, tol=0.0001, Trace=False): > > """ > > Calculate rational approximation of x to within tolerance of tol; > > returns a tuple consisting of numerator and denominator p/q > > Trace=True causes iterated results to be shown > > """ > > a, r, p, q = [], [], [], [] > > Done = False > > n = 0 > > if Trace: print "x:%f tol:%f" % (x, tol) > > while not Done: > > a.append(None) > > r.append(None) > > p.append(None) > > q.append(None) > > if n == 0: r[n] = x > > else: r[n] = 1/(r[n-1]-a[n-1]) > > a[n] = int(r[n]) > > if n == 0: > > p[n] = a[0] > > q[n] = 1 > > elif n ==1: > > p[n] = a[n]*p[n-1] + 1 > > q[n] = a[n] > > else: > > p[n] = a[n]*p[n-1] + p[n-2] > > q[n] = a[n]*q[n-1] + q[n-2] > > if Trace: > > print "n:%d a:%d p:%d q:%d approx:%f" % \ > > (n, a[n], p[n], q[n], float(p[n])/q[n]) > > if abs(float(p[n])/q[n] - x) < tol: > > Done = True > > num = p[n]; denom = q[n] > > n += 1 > > return (num, denom) > > > >##################################################### > > > >Here's a result for pi: > > > > >>> print cf(3.14159265357989,0.0000001, Trace=True) > >x:3.141593 tol:0.000000 > >n:0 a:3 p:3 q:1 approx:3.000000 > >n:1 a:7 p:22 q:7 approx:3.142857 > >n:2 a:15 p:333 q:106 approx:3.141509 > >n:3 a:1 p:355 q:113 approx:3.141593 > >n:4 a:292 p:103993 q:33102 approx:3.141593 > >(103993, 33102) > > > >i.e., the first 5 approximations it came up with were 3/1, 22/7, 333/106, > >355/113 and a whopping 103993/33102. > > > >For the 0.36 example you used earlier: > > > > >>> print cf(0.36, .01, Trace= True) > >x:0.360000 tol:0.010000 > >n:0 a:0 p:0 q:1 approx:0.000000 > >n:1 a:2 p:1 q:2 approx:0.500000 > >n:2 a:1 p:1 q:3 approx:0.333333 > >n:3 a:3 p:4 q:11 approx:0.363636 > >(4, 11) > > >>> > > > >it went right from 1/3 to 4/11 (0.363636), skipping the 3/8 (0.375) option > >from the farey series. > > > >But this continued fraction algorithm is ill-suited to answer the question > >"what's the closest fraction with a denominator < N", because it doesn't > >try to find that, it jumps further ahead with each iteration. > > > >Anyway, I thought you might find it interesting based on our discussion. > >Terry, > >Well, I have to admit I don't understand your code at all. But I see it works. > >I modified one of my functions of frac.py, and came up with > >=============================================== >from __future__ import division >import time, psyco > >psyco.full() > >def d(number): > import decimal > decimal.getcontext().prec = 16 > return decimal.Decimal(str(number)) > >def bestFracForMinimumError(decimal, minimumError): > denom = 0 > smallestError = 10 > count = 0 > while True: > denom += 1 > num = int(round(d(decimal) * d(denom))) > error = abs((((d(num)) / d(denom)) - d(decimal)) / >d(decimal)) * d(100) > if d(error) <= d(smallestError): > count += 1 > smallestError = d(error) > q = d(num)/d(denom) > print "%d/%d = %s has error of %s per cent" % (num, >denom, q, smallestError) > if d(smallestError) <= d(minimumError): > print "count is", count > break > >=====================================================================
I just realized that I had used "<=" in my code where I should have used "<". So after make these changes, the results also changed, of course. See them at http://www.rcblue.com/Python/PartOfReplyToTerryOnTutorList-2.txt Dick >You can see the results of both >bestFracForMinimumError(3.14159265357989, 0.00000002) > >(BTW your pi is a bit off but I used yours, instead of math.pi, which >is 3.1415926535897931 . Also, I needed 0.00000002 in order to produce >your 103993/33102) > >and > >bestFracForMinimumError(.36, .01) > >at <http://www.rcblue.com/Python/PartOfReplyToTerryOnTutorList.txt> _______________________________________________ Tutor maillist - Tutor@python.org http://mail.python.org/mailman/listinfo/tutor
