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

[Terry Carroll]

On Sat, 6 Jan 2007, Dick Moores wrote:

> Well, I have to admit I don't understand your code at all. But I see it
> works.

A continuing fraction is an expression that looks like this:

                1
A +  ----------------------------------
                   1
      B + -----------------------------
                         1
             C + ----------------------
                            1
                  D + -----------------
                         E + etc.

i.e., A + (1 / (B + (1 / (C + (1 / (D + (1 / E + ... )))))))

That is, the denominator of the fraction itself has a term that is a 
fraction; and that fraction itself has a term that is a fraction, etc., 
either to infinity or to some finite point.

The basic idea of using a continued fraction to get a rational 
approximation is something like this.   Say you want to get an 
approximation of 0.096.  Well, a good first approximation is 1/10.  But 
1/10 is 0.010.  That's too big, or, to put another way, the denominator is 
too small.  It should really be 1/(10+something), where "something" is 
some fraction 1/x (if it were larger than 1, then we'd be starting from 
11+ something).

The next best approximation is to figure out that "something" is.   It 
turns out, it's 1/2, i.e. .096 = about 1/(10+(1/2)), or about 0.9624.  
That's still too big, so that 1/2 we added should really be 1/(2+ 
something).

The next "something" we calculate is again 1/2, so it's 
1/(10+(1/(2+(1/2)), or .096154.  Getting closer, but the same thing.

Again, the next something is 1/2, so we get 1/(10+(1/(2+(1/(2+1/2))), 
which is 0.096, exactly.

So, in that butt-ugly ascii I have above, A=10, B = 2, C = 2 and D = 2,
and it turns out not to have been necessary to go any further.

(I admit I cheated above and picked 0.096 as an example, because some 
other examples, like .093 and .098 actually have a "something" that's 
equal to 1).

Now, I actually don't understand all of the math in the article at 
http://mathworld.wolfram.com/ContinuedFraction.html , which I based this 
on; but enough of it is clear that I could basically just take certain 
equations and translate them into Python:


>     a, r, p, q = [], [], [], []

This is a list of terms; "a" is the successive denominator, i.e. compared
to my description above a[0] is A; a[1] is B; a[2] is C, etc.  r[n] is
used to calculate a[n]; it actually could have been dispensed with.  p[n]
and q[n] are the value of the continued fraction after n iterations.

>         a.append(None)
>         r.append(None)
>         p.append(None)
>         q.append(None)

This is just so I can assign to a[n], r[n], etc.; it makes the code more 
readable.

>         if n == 0: r[n] = x
>         else: r[n] = 1/(r[n-1]-a[n-1])

These statements are equations 8 & 9 from the Wolfram page cited above.

>         a[n] = int(r[n])

This is equation 10.

>         if n == 0:
>             p[n] = a[0]
>             q[n] = 1

These are effectively equations 25 & 26 for the special case of n == 0, 
i.e. equation 24.

>         elif n ==1:
>             p[n] = a[n]*p[n-1] + 1
>             q[n] = a[n]

These are equations 25 & 26 for the special case of n == 1, taking into 
account equations 22 to get values for what would have otherwise looking 
for p[-1], i.e. an entry before the first element in the list.

>         else:
>             p[n] = a[n]*p[n-1] + p[n-2]
>             q[n] = a[n]*q[n-1] + q[n-2]

This is the general case for equations 25 & 26.

The rest is just cranking through the iterations.

> (BTW your pi is a bit off but I used yours, instead of math.pi, which 
> is 3.1415926535897931 . 

Oops.  I transposed the "79" and "89"; that's what I get for going from 
memory.

I may add this algorithm to the cookbook.