[Edu-sig] RE: [Tutor] Need more precise digits

[Tim Peters]

[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