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/