[Edu-sig] RE: [Tutor] Need more precise digits
Tim Peters
tim.one@home.com
Sun, 2 Dec 2001 04:57:49 -0500
[Kirby Urner]
> ...
> By the way, another application of the generator feature
> came up for me on another list recently. I was chatting
> with a community college prof about continued fractions,
> which have the standard form:
>
> q0 + 1
> --------
> q1 + 1
> ----
> q2 + ....
>
> Where q0, q1... are positive integers. I had this recursive
> way of doing it which amounted to a right to left evaluation,
> but he showed me a way of going left to right. That was cool,
> because one thing I wanted to see was how some continued
> fractions gradually converge to a key constant. For example,
> the simplest of all continued fractions is just:
>
> 1 + 1
> --------
> 1 + 1
> ----
> 1 + ....
>
> and so on.
Note that the numerators and denominators of the convergents to this are the
Fibonacci numbers; that is, this:
> ...
> 1/1
> 2/1
> 3/2
> 5/3
> 8/5
> 13/8
> 21/13
> ...
should look familiar <wink>.
> The professor had written his algorithm for the TI calculator,
> but I adapted with hardly any modifications to Python, and
> now it looks like this:
>
> from __future__ import generators
>
> def conv(cf):
> """
> Generate partial fractions from partial quotients
> of a continued fraction, with thanks to:
> RWW Taylor
> National Technical Institute for the Deaf
> Rochester Institute of Technology
> Rochester NY 14623
> """
> cfsize = len(cf)
> n = [0 for i in range(cfsize+2)] # all 0s
Simpler as
n = [0] * (cfsize + 2)
but, as shown later, you don't need a list here at all.
> d = n[:] # copy of n
> n[1] = d[0] = 1
> for i in range(cfsize):
> n[i+2] = n[i] + cf[i] * n[i + 1]
> d[i+2] = d[i] + cf[i] * d[i + 1]
> yield str(n[i])+"/"+str(d[i]) # interim result
Consider the last iteration of the loop: the n[i+2]/d[i+2] and
n[i+1]/d[i+1] convergents are computed but never returned (when
i==len(cfsize)-1).
> return
Not needed -- "falling off the end" of a generator is the same as a
"return".
> So if I feed this a list like [1,1,1,1,1,1,1,1,1,1,1,1,1],
Easier written as [1]*13 (see the similar trick with [0] above).
> it'll yield an approximation with each iteration, because
> I wrote it as a generator.
However, because you compute len(cf) at the start, cf can't *itself* be a
generator. It's more flexible (and the code gets simpler!) if you let cf be
any iterable object. For example, then you could feed it this:
def ones():
while 1:
yield 1
That is, an unbounded sequence of ones.
> For example:
>
> >>> genphi = conv([1,1,1,1,1,1,1,1,1,1,1,1,1])
> >>> for fraction in genphi: print fraction
>
> 0/1
> 1/0
> 1/1
> 2/1
> 3/2
> 5/3
> 8/5
> 13/8
> 21/13
> 34/21
> 55/34
> 89/55
> 144/89
>
> Of course a 'print' statement in place of yield would accomplish
> the same thing in this context (that's what a yield is like,
> a print). But I can imagine an application where I'd want
> to return the partial fraction, then maybe ask for more
> precision later on, and conv() would be there, ready to pick
> up right where it left off.
Yes indeed! There are many applications for this, although it's hard to
give an obvious example <wink>.
Here's an alternative that accepts any iterable object (incl. a generator,
if you like) as argument, generates all the convergents (including the last
two), and doesn't use any internal lists:
from __future__ import generators, division
# Generate continued-fraction pairs (num, den) from a sequence of
# partial quotients.
def conv(pqs):
x0, y0 = 0, 1 # zero
x1, y1 = 1, 0 # infinity
yield x0, y0
yield x1, y1
for q in pqs:
x0, y0, x1, y1 = x1, y1, x0 + q*x1, y0 + q*y1
yield x1, y1
import math
x = (1 + math.sqrt(5))/2
for n, d in conv([1] * 50):
approx = n/(d or 1e-300)
print "%d/%d ~= %.17g %.17g" % (n, d, approx, x - approx)
> ...
> This is where we could write a loop that continues until the
> difference between the fraction and the floating point goes
> below some very small value, like 1E-12. Just modify conv()
> to never stop on its own, and keep looping until you hit the
> break condition, which should *not* involve an ==, for the
> reason you mention.
Continued fractions are lovely. One of the things you can prove is that
successive convergents are alternately larger and smaller than the true
value (e.g., 0/1 < phi, 1/0 > phi, 1/1 < phi, 2/1 > phi, <, >, <, >, ...).
Another is that if p/q and r/s are successive convergents, then abs(p/q -
r/s) == 1/(q*s) (e.g., 13/8-21/13 == (13*13-21*8)/(8*13) == 1/(8*13)).
Together, those imply that the true value is within 1/(q*s) of both
convergents. A weaker but more useful relation is that, since the
denominators increase once the sequence gets going, the convergent following
p/q has a denominator at least as large as q, so the "1/(q*s)" is no larger
than 1/q**2. IOW, once you get beyond the trivial convergents at the start,
any particular convergent p/q is within 1/q**2 of the true value. This can
be used to determine a stopping point good to a given level of accuracy even
when you don't know the true value in advance.
One other factoid of use: if p/q is a convergent to a real number x, p/q is
the best rational approximation to x among all rationals with denominator no
greater than q -- although, as usual with continued fraction, that really
needs some weasle words to exempt the trivial 0/1 and 1/0 convergents at the
start.
more-fun-than-apple-pi-ly y'rs - tim