[Tutor] Need help with rewriting script to use Decimal module

[Terry Carroll]

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