back to list

best rational approximation

🔗Carl Lumma <ekin@lumma.org>

5/26/2006 10:36:03 PM

In August of last year I asked about this (where best
is 'lowest Tenney height') and people told me to check
out Euclid's algorithm. Though I know how to use it to
find greatest common divisors, I can't seem to figure
out how to use it to find rational approximations.

Around that time, I posted some code I thought might
do the trick. Here's some output to see if it's right.
The first parameter is the target, the second is the
factor of acceptable error, and the final number is the
proposed answer.

(gear (sqrt 2) 1.4) 3/2
(gear (sqrt 2) 1.07) 3/2
(gear (sqrt 2) 1.06) 7/5
(gear (sqrt 2) 1.02) 7/5
(gear (sqrt 2) 1.01) 17/12
(gear (sqrt 2) 1.001) 41/29
(gear (sqrt 2) 1.0001) 99/70

(gear (sqrt 5) 1.01) 9/4
(gear (sqrt 5) 1.001) 38/17
(gear (sqrt 5) 1.0001) 161/72
(gear (sqrt 5) 1.00001) 521/233
(gear (sqrt 5) 1.000001) 2207/987

-Carl

🔗Gene Ward Smith <genewardsmith@coolgoose.com>

5/26/2006 10:56:56 PM

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:
>
> In August of last year I asked about this (where best
> is 'lowest Tenney height') and people told me to check
> out Euclid's algorithm. Though I know how to use it to
> find greatest common divisors, I can't seem to figure
> out how to use it to find rational approximations.

> (gear (sqrt 5) 1.00001) 521/233
> (gear (sqrt 5) 1.000001) 2207/987

You seem to be getting round off errors, but otherwise are on the
right track.

The Euclidean algorithm stuff has to do with simple continued fractions.
Check out:

http://en.wikipedia.org/wiki/Continued_fraction

Plus Hardy and Wright or Niven and Zuckerman.

🔗Aaron Krister Johnson <aaron@akjmusic.com>

5/27/2006 2:41:27 AM

On Saturday 27 May 2006 12:24 pm, tuning-math@yahoogroups.com wrote:

> In August of last year I asked about this (where best
> is 'lowest Tenney height') and people told me to check
> out Euclid's algorithm. Though I know how to use it to
> find greatest common divisors, I can't seem to figure
> out how to use it to find rational approximations.
>
> Around that time, I posted some code I thought might
> do the trick. Here's some output to see if it's right.
> The first parameter is the target, the second is the
> factor of acceptable error, and the final number is the
> proposed answer.
>
> (gear (sqrt 2) 1.4) 3/2
> (gear (sqrt 2) 1.07) 3/2
> (gear (sqrt 2) 1.06) 7/5
> (gear (sqrt 2) 1.02) 7/5
> (gear (sqrt 2) 1.01) 17/12
> (gear (sqrt 2) 1.001) 41/29
> (gear (sqrt 2) 1.0001) 99/70
>
> (gear (sqrt 5) 1.01) 9/4
> (gear (sqrt 5) 1.001) 38/17
> (gear (sqrt 5) 1.0001) 161/72
> (gear (sqrt 5) 1.00001) 521/233
> (gear (sqrt 5) 1.000001) 2207/987

Carl,

The other, 'brute force' thing you could do is set up the Stern-Brocot tree
and look for approximations to the left and right (greater than or less than)
of your target.

When I have time, I'll post some Python code I've used for this very thing (I
came across this technique last week, actually, while investigating G.
Secor's equal-beating temperaments, and the rationale behind his
approximations)

But I'm sure some kind of algorithm like Gene suggest, using continued
fractions, is more direct (and probably closely related in concept)

Gene, care to post any pseudo-code (or actual code?)

-Aaron.

🔗Carl Lumma <ekin@lumma.org>

5/27/2006 9:50:02 AM

At 02:41 AM 5/27/2006, you wrote:
>On Saturday 27 May 2006 12:24 pm, tuning-math@yahoogroups.com wrote:
>
>> In August of last year I asked about this (where best
>> is 'lowest Tenney height') and people told me to check
>> out Euclid's algorithm. Though I know how to use it to
>> find greatest common divisors, I can't seem to figure
>> out how to use it to find rational approximations.
>>
>> Around that time, I posted some code I thought might
>> do the trick. Here's some output to see if it's right.
>> The first parameter is the target, the second is the
>> factor of acceptable error, and the final number is the
>> proposed answer.
>>
>> (gear (sqrt 2) 1.4) 3/2
>> (gear (sqrt 2) 1.07) 3/2
>> (gear (sqrt 2) 1.06) 7/5
>> (gear (sqrt 2) 1.02) 7/5
>> (gear (sqrt 2) 1.01) 17/12
>> (gear (sqrt 2) 1.001) 41/29
>> (gear (sqrt 2) 1.0001) 99/70
>>
>> (gear (sqrt 5) 1.01) 9/4
>> (gear (sqrt 5) 1.001) 38/17
>> (gear (sqrt 5) 1.0001) 161/72
>> (gear (sqrt 5) 1.00001) 521/233
>> (gear (sqrt 5) 1.000001) 2207/987
>
>Carl,
>
>The other, 'brute force' thing you could do is set up the Stern-Brocot tree
>and look for approximations to the left and right (greater than or less than)
>of your target.

That sounds like what I'm doing!

-Carl

🔗Carl Lumma <ekin@lumma.org>

5/27/2006 10:10:13 AM

>> (gear (sqrt 5) 1.00001) 521/233
>> (gear (sqrt 5) 1.000001) 2207/987
>
>You seem to be getting round off errors, but otherwise are on the
>right track.

I seed the tree with floor(target/range) and ceiling(target*range),
which seems to me the appropriate thing. Here's what I do:

s1 s2
a
b
c
t

S1 and S2 are the seeds, and A is their mediant. Now I take
the mediant between A and the previous generation's number
that's on the other side of the target T from it (in this case
S2, yielding B). Then between B and A yielding C, and let's
say C is within range so I return it. No rounding should be
happening. Here's a trace of mediants used to find the "gear"
from (sqrt 5) within a factor of 1.000001 (which I'm claming
is 2207/987)...

|(middle 2 3)
|5/2
|(middle 2 5/2)
|7/3
|(middle 2 7/3)
|9/4
|(middle 2 9/4)
|11/5
|(middle 11/5 9/4)
|20/9
|(middle 20/9 9/4)
|29/13
|(middle 29/13 9/4)
|38/17
|(middle 38/17 9/4)
|47/21
|(middle 38/17 47/21)
|85/38
|(middle 38/17 85/38)
|123/55
|(middle 38/17 123/55)
|161/72
|(middle 38/17 161/72)
|199/89
|(middle 199/89 161/72)
|360/161
|(middle 360/161 161/72)
|521/233
|(middle 521/233 161/72)
|682/305
|(middle 682/305 161/72)
|843/377
|(middle 682/305 843/377)
|1525/682
|(middle 682/305 1525/682)
|2207/987

>The Euclidean algorithm stuff has to do with simple continued fractions.
>Check out:
>
>http://en.wikipedia.org/wiki/Continued_fraction

I looked at that, but it's damn confusing. However, continued
fractions are apparently paths through the Stern Brocot tree,
so I think if I do it right I'm ok.

>Plus Hardy and Wright or Niven and Zuckerman.

Tx, I'll search for those.

-Carl

🔗Gene Ward Smith <genewardsmith@coolgoose.com>

5/27/2006 10:45:45 AM

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> 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.

🔗Aaron Krister Johnson <aaron@akjmusic.com>

5/27/2006 11:17:52 AM

--- In tuning-math@yahoogroups.com, "Gene Ward Smith"
<genewardsmith@...> wrote:
>
> --- In tuning-math@yahoogroups.com, Carl Lumma <ekin@> 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.
>
Thanks, Gene!

I see the connection between
h[-2] = 0, h[-1] = 1
and
k[-2] = 1, k[-1] = 0
must be the first seeds of the Stern-Brocot tree, yes?
i.e., 0/1 , 1/0, right?

Fascinating!

This is related to the Euclidean GCD algorithm too, of course...

I would love to post some Python code for this, and run it, and you
could confirm that it works.

Carl, I'm still trying to post what you wrote, I think we were
thinking along the same lines, however, I can't understand whether
you included in your idea the understanding that unless one starts
with true Stern-Brocot seed fractions, one could end up with errors
in that the SB-tree should never contain fractions which could be
reduced by a common factor--i.e., the numerators and denominators of
each should be relatively prime. The only way to insure this would
be (apparently) the algorithm Gene gave, or to start the tree with
'0/1 1/0', and make sure that your upper and lower bounds are always
derived as true medians from two adjacent terms from a generation
that derives from that initial seed....(I hope what I just wrote
makes sense to the casual reader!)

-Aaron

🔗Aaron Krister Johnson <aaron@akjmusic.com>

5/27/2006 11:14:24 AM

--- In tuning-math@yahoogroups.com, "Gene Ward Smith"
<genewardsmith@...> wrote:
>
> --- In tuning-math@yahoogroups.com, Carl Lumma <ekin@> 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.

Thanks, Gene!

I see the connection between
h[-2] = 0, h[-1] = 1
and
k[-2] = 1, k[-1] = 0
must be the first seeds of the Stern-Brocot tree, yes?
i.e., 0/1 , 1/0, right?

Fascinating!

This is related to the Euclidean GCD algorithm too, of course...

I would love to post some Python code for this, and run it, and you
could confirm that it works.

Carl, I'm still trying to post what you wrote, I think we were
thinking along the same lines, however, I can't understand whether
you included in your idea the understanding that unless one starts
with true Stern-Brocot seed fractions, one could end up with errors
in that the SB-tree should never contain fractions which could be
reduced by a common factor--i.e., the numerators and denominators of
each should be relatively prime. The only way to insure this would
be (apparently) the algorithm Gene gave, or to start the tree with
'0/1 1/0', and make sure that your upper and lower bounds are always
derived as true medians from two adjacent terms from a generation
that derives from that initial seed....(I hope what I just wrote
makes sense to the casual reader!)

-Aaron

🔗Kees van Prooijen <lists@kees.cc>

5/27/2006 11:54:39 AM

This is also as I used it.

Have a look:

http://www.kees.cc/tuning/interface.html

Kees

--- In tuning-math@yahoogroups.com, "Gene Ward Smith"
<genewardsmith@...> wrote:
>
> --- In tuning-math@yahoogroups.com, Carl Lumma <ekin@> 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.
>

🔗Carl Lumma <ekin@lumma.org>

5/27/2006 2:18:07 PM

>I would love to post some Python code for this, and run it, and you
>could confirm that it works.

Sure, I have Python here.

>Carl, I'm still trying to post what you wrote, I think we were
>thinking along the same lines, however, I can't understand whether
>you included in your idea the understanding that unless one starts
>with true Stern-Brocot seed fractions, one could end up with errors
>in that the SB-tree should never contain fractions which could be
>reduced by a common factor--i.e., the numerators and denominators of
>each should be relatively prime. The only way to insure this would
>be (apparently) the algorithm Gene gave, or to start the tree with
>'0/1 1/0', and make sure that your upper and lower bounds are always
>derived as true medians from two adjacent terms from a generation
>that derives from that initial seed....(I hope what I just wrote
>makes sense to the casual reader!)

Is that true? Brocot's original procedure uses seeds that
are the floor and ceiling of the target.

http://www.americanscientist.org/template/AssetDetail/assetid/20826?print=yes

-Carl

🔗Carl Lumma <ekin@lumma.org>

5/27/2006 3:30:26 PM

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

This looks like what I've been seeing, attributed apparently
to Gauss. Here you've explained it more clearly than I've
seen it elsewhere. I guess I'll code it up. Can you prove
it gives least Tenney Height results (I've only seen that it
gives least denominator results)?

-Carl

🔗Carl Lumma <ekin@lumma.org>

5/29/2006 11:08:30 AM

Since nobody can find a counterexample to the results I've given,
I'll assume my procedure works and there's no need to worry about
Gauss.

-Carl

🔗Gene Ward Smith <genewardsmith@coolgoose.com>

5/29/2006 11:15:51 AM

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:
>
> Since nobody can find a counterexample to the results I've given,
> I'll assume my procedure works and there's no need to worry about
> Gauss.

I thought I pointed out that some of the values were wrong? Though I
think that was round-off error.

🔗Carl Lumma <ekin@lumma.org>

5/29/2006 12:10:24 PM

>> Since nobody can find a counterexample to the results I've given,
>> I'll assume my procedure works and there's no need to worry about
>> Gauss.
>
>I thought I pointed out that some of the values were wrong? Though I
>think that was round-off error.

There was:

"
> (gear (sqrt 5) 1.00001) 521/233
> (gear (sqrt 5) 1.000001) 2207/987

You seem to be getting round off errors, but otherwise are on the
right track.
"

But what are the correct answers?

-Carl

🔗Gene Ward Smith <genewardsmith@coolgoose.com>

5/29/2006 1:40:39 PM

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

> But what are the correct answers?

For sqrt(5), a particularly simple example BTW, it's 2, 9/4, 38/17,
161/72, 682/305, 2889/1292, 12238/5473. If N[i]/D[i] is the nth term,
then the n+1 term is (4*N[i]+N[i-1])/(4*D[i]+D[i-1]), making them easy
to compute. Another fun fact is that if you substitute N[i]/D[i] into
(x+5/x)/2, you get N[2i]/D[2i].

🔗Carl Lumma <ekin@lumma.org>

5/29/2006 1:57:57 PM

>>> (gear (sqrt 5) 1.00001) 521/233
>>> (gear (sqrt 5) 1.000001) 2207/987
>>
>>You seem to be getting round off errors, but otherwise are on the
>>right track.
>
>But what are the correct answers?
>
>For sqrt(5), a particularly simple example BTW, it's 2, 9/4, 38/17,
>161/72, 682/305, 2889/1292, 12238/5473.

Ok, so the first of these that's within a factor of 1.00001
(north or south) of sqrt(5) is 682/305...

> (/ (sqrt 5) 682/305)
1.0000010749815775

161/72 just misses...

> (/ 161/72 (sqrt 5))
1.0000192899374059

And we have

> (tenneyheight 682/305)
208010

Meanwhile, my answer, 521/233, is within the range

> (/ (sqrt 5) 521/233)
1.0000073680565278

and has a *lower tenney height*...

> (tenneyheight 521/233)
121393

For the tighter 1.000001 range, it looks like we have to go to
2889/1292...

> (/ 2889/1292 (sqrt 5))
1.0000000599066396

And again, my answer...

> (/ 2207/987 (sqrt 5))
1.000000410606289

has lower tenney height and denominator...

> (tenneyheight 2889/1292)
3732588

> (tenneyheight 2207/987)
2178309

-Carl

🔗Gene Ward Smith <genewardsmith@coolgoose.com>

5/29/2006 3:10:36 PM

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

> Ok, so the first of these that's within a factor of 1.00001
> (north or south) of sqrt(5) is 682/305...
>
> > (/ (sqrt 5) 682/305)
> 1.0000010749815775

D^2 * |sqrt(5)-N/D| = 0.2236

> 161/72 just misses...
>
> > (/ 161/72 (sqrt 5))
> 1.0000192899374059

D^2 * |sqrt(5)-N/D| = 0.2236

> And we have
>
> > (tenneyheight 682/305)
> 208010

D^2 * |sqrt(5)-N/D| = 0.2236

> Meanwhile, my answer, 521/233, is within the range
>
> > (/ (sqrt 5) 521/233)
> 1.0000073680565278
>
> and has a *lower tenney height*...

D^2 * |sqrt(5)-N/D| = 0.8944 (not as good!)

> For the tighter 1.000001 range, it looks like we have to go to
> 2889/1292...
>
> > (/ 2889/1292 (sqrt 5))
> 1.0000000599066396

D^2 * |sqrt(5)-N/D| = 0.2236

> And again, my answer...
>
> > (/ 2207/987 (sqrt 5))
> 1.000000410606289

D^2 * |sqrt(5)-N/D| = 0.8944

Your approximations are semiconvergents, but not convergents.

🔗Carl Lumma <ekin@lumma.org>

5/29/2006 6:16:38 PM

>> Ok, so the first of these that's within a factor of 1.00001
>> (north or south) of sqrt(5) is 682/305...
>>
>> > (/ (sqrt 5) 682/305)
>> 1.0000010749815775
>
>D^2 * |sqrt(5)-N/D| = 0.2236

What's the significance of D^2 ??

>Your approximations are semiconvergents, but not convergents.

Who cares? I just want the least Tenney height within the
range.

-Carl

🔗Gene Ward Smith <genewardsmith@coolgoose.com>

5/29/2006 6:45:56 PM

--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:
>
> >> Ok, so the first of these that's within a factor of 1.00001
> >> (north or south) of sqrt(5) is 682/305...
> >>
> >> > (/ (sqrt 5) 682/305)
> >> 1.0000010749815775
> >
> >D^2 * |sqrt(5)-N/D| = 0.2236
>
> What's the significance of D^2 ??

I'm just showing that relative to the size of the numbers involved,
your approximations are not the best possible. You could relate this
to Tenney height by multiplying the error by the numerator times the
denominator. If e is the error and T is the Tenney height, we get to
close approximation that Te = constant1 for the best values, and
constant2 for the second best.

🔗Carl Lumma <ekin@lumma.org>

5/29/2006 8:06:07 PM

At 06:45 PM 5/29/2006, you wrote:
>--- In tuning-math@yahoogroups.com, Carl Lumma <ekin@...> wrote:
>>
>> >> Ok, so the first of these that's within a factor of 1.00001
>> >> (north or south) of sqrt(5) is 682/305...
>> >>
>> >> > (/ (sqrt 5) 682/305)
>> >> 1.0000010749815775
>> >
>> >D^2 * |sqrt(5)-N/D| = 0.2236
>>
>> What's the significance of D^2 ??
>
>I'm just showing that relative to the size of the numbers involved,
>your approximations are not the best possible.

I just want the simplest number in the range, nothing more.

-Carl