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

Kirby Urner urnerk@qwest.net
Sun, 02 Dec 2001 19:11:59 -0800


>
>more-fun-than-apple-pi-ly y'rs  - tim

Excellent stuff on continued fractions from our
resident math-through-Python (and Python-thru-math)
guru Tim Peters.

Going the other direction, you might want to take some
ordinary fraction p/q, and convert it into a continued
fraction, expressed as a list of partial quotients
[q0,q1,q2...].

Here's an algorithm for doing that from my algebra.py,
based on stuff I'm pretty sure I learned from another
good math book: 'Number' by Midhat Gazale.

  def cfract(a,b):
      """
      Return partial quotients of a regular
      continued fraction equivalent to a/b"""
      rcf = []
      while b<>0:
         p = a//b
         rcf.append(p)
         b, a = a - p*b, b
      return rcf

It's sort of a modification of the EEA (Euclid's Extended
Algorithm).

So suppose we want the continued fraction expression of phi,
starting from x = (1 + sqrt(5)/2.  You can just put the
floating point decimal over a big power of 10, and run
it through:

  >>> import algebra, math
  >>> phi = (1 + math.sqrt(5))/2
  >>> phi
  1.6180339887498949
  >>> 16180339887498949/10000000000000000.
  1.6180339887498949
  >>> algebra.cfract(16180339887498949,10000000000000000)
  [1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 10L, 8L, 5L, 37L, 1L,
  10L, 1L, 2L, 3L, 3L, 1L, 1L, 1L, 1L, 2L]

OK, so it goes off the rails there after awhile, or
actually, it doesn't, because this is an exact translation
of the above fraction, which isn't phi, but an approximation
of phi.

Another interesting true fact about continued fractions, is
the square root of any natural number will generate a
*repeating pattern* of partial quotients.  Taking this as
a given, we can use the same trick as above to get a handle
on what this pattern might be:

  >>> math.sqrt(3)
  1.7320508075688772
  >>> algebra.cfract(17320508075688772,10000000000000000)[:15]
  [1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L]

I truncated the list to 15 elements, knowing it'll diverge
from the pattern later, given I'm feeding it an approximation.
This is enough to give me the idea.

OK, so it looks like I've got a repeating 1,2 pattern, with
a 1 up front.  I could yield that with a generator, as we've
seen:

   def root3pqs():
       yield 1
       flipper = 1
       while 1:
          if flipper:
              yield 1
              flipper = 0
          else:
              yield 2
              flipper = 1

(that could be made more streamlined I'm sure).  Testing:

   >>> genpqs = root3pqs()
   >>> [genpqs.next() for i in range(10)]
   [1, 1, 2, 1, 2, 1, 2, 1, 2, 1]

Lookin' good.

So now I can use Tim's cf() generator with this generator as
input, to get successively more accurate fractional
representations of the square root of 3:

  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

  >>> genpqs = root3pqs()
  >>> root3gen = conv(genpqs)
  >>> root3gen.next()
  (0, 1)
  >>> root3gen.next()
  (1, 0)
  >>> root3gen.next()
  (1, 1)
  >>> root3gen.next()
  (2, 1)
  >>> root3gen.next()
  (5, 3)
  >>> root3gen.next()
  (7, 4)
  >>> for i in range(30):      # skip ahead in a hurry
         val = root3gen.next()

  >>> val
  (2642885282L, 1525870529)
  >>> 2642885282L/1525870529.
  1.7320508075688772
  >>> math.sqrt(3)
  1.7320508075688772

Our fraction is already quite accurate.

BTW, I second Tim's recommendation of 'Concrete Mathematics'
-- bought it awhile ago on his recommendation and never
regretted it, even though I'm not a CS major (philosophy
over here).

Kirby