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. (Yeah, I know, I didn't choose the formats well on those print statements.) _______________________________________________ Tutor maillist - Tutor@python.org http://mail.python.org/mailman/listinfo/tutor
