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.
_______________________________________________
Tutor maillist - Tutor@python.org
http://mail.python.org/mailman/listinfo/tutor
