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.

>Here's an example of a shell session where I approximated the sqrt(2)
>with rationals:
//
>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
...

Here are the mediants I find on the way to an answer that
is within a factor of 1.0000001 of sqrt(2)...

|3/2
|4/3
|7/5
|10/7
|17/12
|24/17
|41/29
|58/41
|99/70
|140/99
|239/169
|338/239
|577/408
|816/577
|1393/985
|1970/1393
|3363/2378

Your results are every-other number of mine. As I've showed,
this lets me find simpler ratios within a range of the target.
It's this range bit that everyone here seems to be ignoring,
but which is a central part of music theory applications from
harmonic entropy to beat-synched scales. Whether I'm missing
any ratios is a question I could use help answering.

-Carl

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:

> It's this range bit that everyone here seems to be ignoring,
> but which is a central part of music theory applications from
> harmonic entropy to beat-synched scales. Whether I'm missing
> any ratios is a question I could use help answering.

And the question is?

Aaron Krister Johnson wrote:

> 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

By the by, I've decided to rewrite this to use generators and not worry about keeping track of the numerators:

from __future__ import generators
from __future__ import division

empirical_accuracy_limit = 24

def ra(x):
"""nice output for the rational approximation"""
for pair in convergent(x):
print "%i/%i" % pair

def convergent(x):
"""'rational approximation' or convergent"""
prev, curr = 1, 0

for num in contfrac(x):
prev, curr = curr, num*curr + prev
yield nint(x*curr), curr

def contfrac(a):
"""Continued fraction expansion"""
b=1
for count in range(empirical_accuracy_limit):
if b==0:
return
yield int(a/b)
a, b = b, a%b

def nint(x):
"""nearest integer"""
return int(round(x))

It works in 2.2 and 2.3. You shouldn't need the future imports in 2.4 because it has generators as standard and the division shouldn't matter (see below)

> 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

Here's my equivalent

>>> ra(math.sqrt(2))
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

I think it's the same

> > >>>>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.

You can also use the future division, so that / always returns floating point.

>>> from __future__ import division
>>> 209064253/147830751
1.4142135623730951

In your code it's better to do float(b) than (b+0.0) because it's clearer what your intention is. Better would be to initialize b as a float. But even then I don't think you need to worry. If a is an integer you'll only get 1 result with b as 1. If a is a float then everything gets promoted to floats automatically.

Also note the whole expression

floor(a/(b+0.0))

which returns a float could be simplified to

a//b

with future division. But I prefer to return an integer.

Graham

>> It's this range bit that everyone here seems to be ignoring,
>> but which is a central part of music theory applications from
>> harmonic entropy to beat-synched scales. Whether I'm missing
>> any ratios is a question I could use help answering.
>
>And the question is?

I tried to be very clear about the problem I'm trying to solve --
that of finding ratios of least tenney height within an arbitrary
range of an arbitrary target. I've shown convergents don't do the
job. Does the procedure I've described?

-Carl

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:

> I tried to be very clear about the problem I'm trying to solve --
> that of finding ratios of least tenney height within an arbitrary
> range of an arbitrary target. I've shown convergents don't do the
> job. Does the procedure I've described?

What would do the job is pretty simple--look at all the
semiconvergents up to your range cutoff.

>> I tried to be very clear about the problem I'm trying to solve --
>> that of finding ratios of least tenney height within an arbitrary
>> range of an arbitrary target. I've shown convergents don't do the
>> job. Does the procedure I've described?
>
>What would do the job is pretty simple--look at all the
>semiconvergents up to your range cutoff.

Thanks! I think that's what I'm doing. I'll have to check.

-Carl

>> I tried to be very clear about the problem I'm trying to solve --
>> that of finding ratios of least tenney height within an arbitrary
>> range of an arbitrary target. I've shown convergents don't do the
>> job. Does the procedure I've described?
>
>What would do the job is pretty simple--look at all the
>semiconvergents up to your range cutoff.

Ok, so here's what Wikipedia says about them:

>If h_(n-1)/k_(n-1) and h_n/k_n are successive convergents,
>then any fraction of the form
>
>h_(n-1) + ah_n
>--------------
>k_(n-1) + ak_n
>
>where 'a' is a nonnegative integer and the numerators and
>denominators are between the n and n + 1 terms inclusive are called
>semiconvergents, secondary convergents, or intermediate fractions.
>Often the term is taken to mean that being a semiconvergent excludes
>the possibility of being a convergent, rather than that a convergent
>is a kind of semiconvergent.
>The semiconvergents to the continued fraction expansion of a real
>number x include all the rational approximations which are better
>than any approximation with a smaller denominator. Another useful
>property is that consecutive semiconvergents a/b and c/d are such
>that ad - bc = +/-1.

Did you mean the term to exclude or include convergents?
My results...

>Here are the mediants I find on the way to an answer that
>is within a factor of 1.0000001 of sqrt(2)...
>
>|3/2
>|4/3
>|7/5
>|10/7
>|17/12
>|24/17
>|41/29
>|58/41
>|99/70
>|140/99
>|239/169
>|338/239
>|577/408
>|816/577
>|1393/985
>|1970/1393
>|3363/2378

...appear to meet the ad - bc = +/-1 condition.

-Carl

>Did you mean the term to exclude or include convergents?
>My results...
>
>>Here are the mediants I find on the way to an answer that
>>is within a factor of 1.0000001 of sqrt(2)...
>>
>>|3/2
>>|4/3
>>|7/5
>>|10/7
>>|17/12
>>|24/17
>>|41/29
>>|58/41
>>|99/70
>>|140/99
>>|239/169
>>|338/239
>>|577/408
>>|816/577
>>|1393/985
>>|1970/1393
>>|3363/2378
>
>...appear to meet the ad - bc = +/-1 condition.
>
>-Carl

Scala will find both convergents and semiconvergents with
its "convergents /semi" command, and even tell you which is which.
My Scheme code seems to match Scala's output here perfectly,
returning both convergents and semiconvergents. That's good.

What's weird is, all the stuff on convergents I read on the
web, hardly anyone mentioned semiconvergents, and nobody mentioned
my algorithm, which is dead simple both conceptually and
computationally.

I still don't have a proof that this catches all ratios with least
n*d, but it seems like it does.

-Carl

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:

> I still don't have a proof that this catches all ratios with least
> n*d, but it seems like it does.

Semiconvergents do that.

>> I still don't have a proof that this catches all ratios with least
>> n*d, but it seems like it does.
>
>Semiconvergents do that.

Is the proof known? Googling for "Hardy Wright semiconvergent"
and "Niven Zuckerman semiconvergent" isn't giving me much luck.

-Carl

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:
>
> >> I still don't have a proof that this catches all ratios with least
> >> n*d, but it seems like it does.
> >
> >Semiconvergents do that.
>
> Is the proof known? Googling for "Hardy Wright semiconvergent"
> and "Niven Zuckerman semiconvergent" isn't giving me much luck.

Niven and Zuckerman stick it in the exercises as a starred (difficult)
exercise, which is sort of annoying. Maybe I need to put the basic
information on semiconvergents on Wikipedia somewhere.

NZ define (in their exercises) a "good approximation" a/b to x as
being one where |bx-a| is minimal over all u/v with v <= b. Then one
exercise is to show an approximation is "good" iff a/b is a
convergent. Then they define a "fair approximation" as one where
|x-a/b| is minimal among all u/v with v <= b, and then ask you to show
that every "fair approximation" is a convergent or a semiconvergent.
Niven and Zukerman is where I learned these facts long, long ago and
not far away in Berkeley as an undergraduate, and I don't know another
source offhand.

Not all semiconvergents are "fair approximations"; 17-et has a fifth
3.93 cents sharp, sharper than 12 is flat. On the other hand, 29-et is
1.49 cents sharper, which beats out 12 and makes it "fair". It isn't
"good" however, because 12-et is 1.955 relative cents flat, and 29-et
is 3.609 relative cents sharp (to get relative cents, multiply the
error in cents by n/12.) The next "good" one is 41-et, which has a
fifth 1.654 relative cents sharp.

The "fair" list for fifths goes 1, 2, 3, 5, 7, 12, 29, 41, 53, 200,
253, 306, 359, 665 ...

--- In tuning-math@yahoogroups.com, "Gene Ward Smith"
<genewardsmith@...> wrote:

> Maybe I need to put the basic
> information on semiconvergents on Wikipedia somewhere.

Actually I did put a section on semiconvergents in there almost a year
ago to the day, which is still there; there is now a section following
it on best rational approximations which I edited to say that
convergents and semiconvergents give you all of the best rational
approximations, which wasn't stated explicitly. So Wikipedia seems to
be doing OK.