As promised, some Python code for continued fraction convergents

> From: "Gene Ward Smith" genewardsmith@coolgoose.com
> Date: Sat May 27, 2006 10:46am(PDT)
> Subject: Re: best rational approximation
>
> Gene wrote:
>
> If you have some x, call it x0. Then set n0 = floor(x0), and
> x1 = 1/(x0-n0). Iterate this, n_i = floor(x_i), and
> x_(i+1) = 1/(x_i-n_i). You then get a sequence of positive integers
> <n0, n1, n2, ...> which is the continued fraction for x; for every
> irrational x, there is a *unique* continued fraction, which is pretty
> neat.
>
> Now set h[-2] = 0, h[-1] = 1, h[i] = n_i h[i-1] + h[i-2], and similarly
> k[-2] = 1, k[-1] = 0, k[i] = n_i k[i-1] + k[i-2]. Then Q[i] = h[i]/k[i]
> is the ith continued fraction approximate for x, starting from i=0.

I'm posting this to the normal tuning groups as well as tuning-math, because I
think it should be dymystified, and it's quite fun and useful!

Hopefully, you will understand what Gene wrote above in a very concrete way by
the time you are done reading this...

This is cool...I figured out that on a calculator with a '1/x' button, to get
the <n0,n1,n2,....> sequence for the simple continued fraction of x, you:

1. Write down the integer part of x (floor(x))
2. Take the integer part of x away from x (x-floor(x))
3. you now have a number between 0 and 1. Take its' inverse
(i.e. '1/x' on the calculator)
4. Goto step '1' and repeat as long as desired...

The caveat is that a calculator, and most computer situations actually, have
limited decimal point precision. So for instance, the continued fraction
expansion for sqrt(2) should be <1,2,2,2,2,2,2,2,2,2,......>, however, once
the precision runs out on a calculator, you will see that chain of '2's gets
disrupted......

Here are some meantone fifths that I figured out by hand with a calculator
using Gene's algorithm; some Python code do do this 'automagically' follows
below. To understand how I got the convergents' numerators and denominators,
take the current continued fraction expansion number, multiply it by the
previous numerator (or denominator) and add the previous-previous numerator
(or denominator).....for example, below, where the top row says '39', we do
(39*3)+1=118
(39*2)+1=79
..which gives us the next convergent, 118/79.

1/3-comma fifth:

1 2 39 1 ;;;; the cont. fraction expans of 1.49380158....

0 1 1 3 118 121 ;;;; the convergents (numerators)
- - - - --- ---
1 0 1 2 79 81 ;;;; the convergents (denoominators)

2/7-comma fifth:

1 2 46 1

0 1 1 3 139 142
- - - - --- ---
1 0 1 2 93 95

1/4-comma fifth:

1 2 53 4

0 1 1 3 160 643
- - - - --- ---
1 0 1 2 107 430

1/5-comma fifth:

1 2 66 1

0 1 1 3 199 202
- - - - --- ---
1 0 1 2 133 135

1/6-comma fifth:

1 2 80 12

0 1 1 3 241 2895
- - - - --- ----
1 0 1 2 161 1934

1/7-comma fifth:

1 2 93 2

0 1 1 3 280 563
- - - - --- ---
1 0 1 2 187 376

1/8-comma fifth:

1 2 106 1

0 1 1 3 319 322
- - - - --- ---
1 0 1 2 213 215

etc, etc........

here's some python code, to be put in a library.
I call my own math functions library 'math_stuff.py'

from math import sqrt, log, floor

def contfrac(a): ### continued fraction expansion
terms=[]
count=0
b=1
while ((b != 0) and (count < 24)): ### 24 times, emprical accuracy
### limit
terms.append(floor(a/(b+0.0)))
a,b = b, a % b
count = count + 1
return terms

def ra(x): ### 'rational approximation', or convergent
numerators=[0,1]
denominators=[1,0]
expansion=contfrac(x) ### call the contfrac function
for num in expansion: ### [-1] and [-2] index 'previous'
### and 'previous-previous'
numerators.append((num*numerators[-1])+numerators[-2])
denominators.append((num*denominators[-1])+denominators[-2])
for index in range(len(numerators)):
print "%i/%i" % (numerators[index], denominators[index])
print

Here's an example of a shell session where I approximated the sqrt(2) with
rationals:

akj@aaron:~$ python
Python 2.4.1 (#1, Apr 10 2005, 22:30:36)
[GCC 3.3.5] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from math_stuff import *
>>> ra(sqrt(2))
0/1
1/0
1/1
3/2
7/5
17/12
41/29
99/70
239/169
577/408
1393/985
3363/2378
8119/5741
19601/13860
47321/33461
114243/80782
275807/195025
665857/470832
1607521/1136689
3880899/2744210
9369319/6625109
22619537/15994428
54608393/38613965
77227930/54608393
131836323/93222358
209064253/147830751

>>> 209064253/147830751.0
1.4142135623730951
>>> sqrt(2)
1.4142135623730951 #### they match!!!
>>>

Notice the use of the '.0' in Python, to tell the interpreter you want a
floating point answer.

Hope this was enlightening!

Best,
Aaron.

--- In tuning@yahoogroups.com, Aaron Krister Johnson <aaron@...> wrote:

> This is cool...I figured out that on a calculator with a '1/x'
button, to get
> the <n0,n1,n2,....> sequence for the simple continued fraction of x,
you:
>
> 1. Write down the integer part of x (floor(x))
> 2. Take the integer part of x away from x (x-floor(x))
> 3. you now have a number between 0 and 1. Take its' inverse
> (i.e. '1/x' on the calculator)
> 4. Goto step '1' and repeat as long as desired...

It's even neater if your calculator has a "frac" key, as some do. That
takes the fractional part, x - floor(x). Then you simply punch the
frac key followed by the 1/x key to get the next value.

It's also useful to know that if you simply want the denominators, as
is sometimes the case, then you don't need the numerators, since
calculating numerators and denominators is a separate calculation.
Also, from just the denominators, you can get the whole sequence as an
infinite series, which also tells you about the closeness of the
approximation.

If I get just the denominators for log2(5), I get a seqence of equal
temperaments for 1/4 comma meantone: 1, 2, 5, 7, 12, 19, 31, 174,
205... Now I can cook up an infinite series from these denominators:

1 - 1/(1*2) + 1/(2*5) - 1/(5*7) + 1/(7*12) - 1/(12*19) + 1/(19*31) ...

The partial sums of this series are simply the convergents of the
continued fraction, which since it is an alternating series tells me
that using 19 as an approximation for 1/4-comma will give an error
less than the next term, 1/(19*31); so the error for 19 is less than
1200/(19*31) = 2.04 cents. In fact the actual error is 1.84 cents. I
also know 19 is flat of the correct value, 31 sharp, and so forth.

--- In tuning@yahoogroups.com, "Gene Ward Smith" <genewardsmith@...>
wrote:
>
> --- In tuning@yahoogroups.com, Aaron Krister Johnson <aaron@>
wrote:
>
> > This is cool...I figured out that on a calculator with a '1/x'
> button, to get
> > the <n0,n1,n2,....> sequence for the simple continued fraction
of x,
> you:
> >
> > 1. Write down the integer part of x (floor(x))
> > 2. Take the integer part of x away from x (x-floor(x))
> > 3. you now have a number between 0 and 1. Take its' inverse
> > (i.e. '1/x' on the calculator)
> > 4. Goto step '1' and repeat as long as desired...
>
> It's even neater if your calculator has a "frac" key, as some do.
That
> takes the fractional part, x - floor(x). Then you simply punch the
> frac key followed by the 1/x key to get the next value.
>
> It's also useful to know that if you simply want the denominators,
as
> is sometimes the case, then you don't need the numerators, since
> calculating numerators and denominators is a separate calculation.
> Also, from just the denominators, you can get the whole sequence
as an
> infinite series, which also tells you about the closeness of the
> approximation.
>
> If I get just the denominators for log2(5), I get a seqence of
equal
> temperaments for 1/4 comma meantone: 1, 2, 5, 7, 12, 19, 31, 174,
> 205... Now I can cook up an infinite series from these
denominators:
>
> 1 - 1/(1*2) + 1/(2*5) - 1/(5*7) + 1/(7*12) - 1/(12*19) +
1/(19*31) ...
>
> The partial sums of this series are simply the convergents of the
> continued fraction, which since it is an alternating series tells
me
> that using 19 as an approximation for 1/4-comma will give an error
> less than the next term, 1/(19*31); so the error for 19 is less
than
> 1200/(19*31) = 2.04 cents. In fact the actual error is 1.84 cents.
I
> also know 19 is flat of the correct value, 31 sharp, and so forth.

Very cool!

-Aaron