best rational approximation

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

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

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.

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

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

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

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

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

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

>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

>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

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

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

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

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

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

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

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

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

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