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 <kovarik@mcmail.cis.mcmaster.ca> 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<d is
the maximum number that will not produce a denominator that is too
large.  These two approximations sandwich the input and are the best
possible with the given denominator bound.  A similar routine with
a different computation of x would give the prescribed-error version
of the best approximation.

For instance with input pi and a denominator bound of 64, the program
below produces as output the fractions 22/7 and 201/64, which are the
closest fractions on either side of pi with such small denominators.
The sequence of approximations on the 201/64 side continues 223/71,
245/78, etc., until one reaches the next convergent.
-- 
David Eppstein		UC Irvine Dept. of Information & Computer Science
eppstein@ics.uci.edu	http://www.ics.uci.edu/~eppstein/

            ------------

/*
** find rational approximation to given real number
** David Eppstein / UC Irvine / 8 Aug 1993
**
** usage: a.out r d
**   r is real number to approx
**   d is the maximum denominator allowed
**
** based on the theory of continued fractions
** if x = a1 + 1/(a2 + 1/(a3 + 1/(a4 + ...)))
** then best approximation is found by truncating this series
** (with some adjustments in the last term).
**
** Note the fraction can be recovered as the first column of the matrix
**  ( a1 1 ) ( a2 1 ) ( a3 1 ) ...
**  ( 1  0 ) ( 1  0 ) ( 1  0 )
** Instead of keeping the sequence of continued fraction terms,
** we just keep the last partial product of these matrices.
*/

#include <stdio.h>

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/