From: eppstein@wormwood.ics.uci.edu (David Eppstein)
Newsgroups: sci.math
Subject: Re: given decimal expansion, find fraction ?
Date: 22 Mar 1995 10:27:19 -0800
In <3kkrnr$klv@apple.com> rjohnson@apple.com (Robert Johnson) writes:
< In article <3j6vld$haf@mcmail.CIS.McMaster.CA>,
< Zdislav V. Kovarik wrote:
< >(And the suspicion that the convergents give approximations with the
< >smallest denominator, subject to prescribed error, is a provable
< >statement.)
< That's interesting, since it's false. Consider everyone's favorite
< transcentdental, pi, and a prescribed error of .001: 22/7 doesn't do
< it, but 333/106 does. However, 201/64 is within .001 of pi and it has
< a smaller denominator.
It is not true (as you point out) that the convergents are always the
best approximations (either with prescribed error or equivalently with a
prescribed bound on the denominator). However the best approximation
can always be found by a small perturbation of the convergents, much
more quickly than your suggestion of using the Stern-Brocot tree.
I include below a short program I wrote a while back for computing best
rational approximations. Similar routines are built into Mathematica
and other symbolic algebra systems. Basically it computes the continued
fraction (a,b,c,d,...), truncates it just before the denominators would
grow too large (giving say (a,b,c)), and then either takes the
convergent of this sequence or of the sequence (a,b,c,x) where x
main(ac, av)
int ac;
char ** av;
{
double atof();
int atoi();
void exit();
long m[2][2];
double x, startx;
long maxden;
long ai;
/* read command line arguments */
if (ac != 3) {
fprintf(stderr, "usage: %s r d\n");
exit(1);
}
startx = x = atof(av[1]);
maxden = atoi(av[2]);
/* initialize matrix */
m[0][0] = m[1][1] = 1;
m[0][1] = m[1][0] = 0;
/* loop finding terms until denom gets too big */
while (m[1][0] * (ai = (long) x) + m[1][1] <= maxden) {
long t;
t = m[0][0] * ai + m[0][1];
m[0][1] = m[0][0];
m[0][0] = t;
t = m[1][0] * ai + m[1][1];
m[1][1] = m[1][0];
m[1][0] = t;
x = 1/(x - (double) ai);
}
/* now remaining x is between 0 and 1/ai */
/* approx as either 0 or 1/m where m is max that will fit in maxden */
/* first try zero */
printf("%ld/%ld, error = %e\n", m[0][0], m[1][0],
startx - ((double) m[0][0] / (double) m[1][0]));
/* now try other possibility */
ai = (maxden - m[1][1]) / m[1][0];
m[0][0] = m[0][0] * ai + m[0][1];
m[1][0] = m[1][0] * ai + m[1][1];
printf("%ld/%ld, error = %e\n", m[0][0], m[1][0],
startx - ((double) m[0][0] / (double) m[1][0]));
}
--
David Eppstein UC Irvine Dept. of Information & Computer Science
eppstein@ics.uci.edu http://www.ics.uci.edu/~eppstein/