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

[Dick Moores]

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>