#!/home/johnh/BIN/perl5 -w

#
# rational_approximations.pl
# $Id$
#
# Copyright (C) 1998-1999 by John Heidemann
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License,
# version 2, as published by the Free Software Foundation.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
#



sub usage {
    print STDERR <<END;
usage: $0 [-v] [denominator-limit]

Outputs a list of rational approximations to pi.
By default, data is shown only when there is a new best approximation.
The verbose flag (-v) shows data for each new denominator.
The debug flag (-d) shows additional information.

The optional denominator-limit specifics how far to compute.

The algorithm is to start with a rational approximation
(initially 3/1) and refine that approximation by testing each
successive denominator.

We make use of the fact that is we know a rational number a/b
which we want to use to approximate R,
then we know that the closest approximation with the next higher
denominator is c/(b+1)
which is one of (c'-1)/(b+1), c'/(b+1), (c'+1)/(b+1)
where c' = a + int(a/b).

Output format is a table of the format
	#h a     b        r                  e                     new
	3        1        3.0000000000000000 1.4159265358979312e-01 *
	13       4        3.2500000000000000 1.0840734641020688e-01 *
	16       5        3.2000000000000002 5.8407346410207062e-02 *
	19       6        3.1666666666666665 2.5074013076873403e-02 *
	22       7        3.1428571428571428 1.2644892673496777e-03 *

where each column is a different value.
The first two (a, b) are the numerator and denominator
of the rational approximation.
The third (r) is the decimal expansion of a/b.
The fourth (e) is the error, i.e. abs(pi - a/b).
The fifth (new) is a flag indicating when the approximation is a new best.

The algorithm is not written using numerically sound techniques, so
its accuracy is limited (although this has not yet been a limit to the
program's amusement value.)  The problem is that it depends on the
accuracy of subtracting two floating point numbers.  Since floating
point numbers are of finite precision (unlike the mathematical concept
of real numbers that they approximate), this means that you cannot
trust this program to approximate pi closer than the number of digits
(minus one) present in your computer's implementation of floating
point.  For IEEE floating point present in all modern computers gives
you at least 14 decimal digits, although I should look up (or get the
energy to compute) exactly how accurate it is.  Eventually the program
should be re-written using a bignum package to get around this
problem.

Another comment about the algorithm:  currently we iterate through all
denominators.  This seems easy and works well since computers are
fast.  It might be better instead to use the error value reported to 
jump more quickly to the next plausible candidate.  (I.e., 22/7ths is
accurate to ~1.3e-3, and numbers expressed in units of a/7 have a
resolution of ~1.4e-1.  From these two we should be able to figure out
that we can't possibly do better for a while.  I haven't worked out
the exact relationship yet.)

END
    #'
    exit 1;
}

use strict;

use Getopt::Long;
&usage if ($#ARGV >= 0 && $ARGV[0] eq '-?');
my(%opts);
&GetOptions(\%opts, qw(d v));
# &usage if ($#ARGV < 0);

my($verbose) = defined($opts{'v'});
my($debug) = defined($opts{'d'});
my($limit) = ($#ARGV >= 0) ? $ARGV[0] : (2<<31);

#
# Start with pi to 40 digits.
# Verified against http://cad.ucla.edu/repository/useful/PI.txt,
# and we know that everything on the web is true.
#
# Note: current floating point systems will discard many of these spare digits.
#
my($R) = 3.1415926535897932384626433832795028841971;
my($a, $e);
my($last_a, $last_e);
$last_a = $a = 3;
$last_e = $R - $last_a;

print "#h a     b        r                  e                     new\n";

my($best_a, $best_b, $best_e);

my($i) = 1;
for (; $i <= $limit;) {
    #
    # show current iteration
    #
    my($r) = $a / ($i + 0.0);
    $e = abs($R - $r);
    my($value) = $e == $last_e ? "same" : ($e < $last_e ? "better" : "worse");

    my($show, $new);
    if (!defined($best_e) || $e < $best_e) {
	$best_a = $a;
	$best_b = $i;
	$best_e = $e;
	$show = 1;
	$new = '*';
    } else {
	$new = '.';
    };

    # print "$a[$i]/$i = $r ($e[$i], $value)\n";n
    printf ("%-8s %-8s %1.16f %1.16e %s\n", $a, $i, $r, $e, $new) if ($show || $verbose);

    #
    # find next iteration
    #
    my($j) = $i+1;
    
    # a/b =  c/b+1
    # 
    # c = a(b+1)/b = ab/b + a/b = a + a/b
    # 
    # ie
    # 
    # 3/1 = 3
    # 6/2 = 3
    #  /3
    my($c) = $a + int($a/$i);
    my($best_d, $best_d_error);
    my($d);
    foreach $d (0, -1, 1) {
	my($trial) = $c + $d;
	my($result) = $trial/($j+0.0);
	my($error) = $R - $result; $error = -$error if ($error < 0);
	my($best) = ' ';
	if (!defined($best_d_error) || $error < $best_d_error) {
	    $best_d = $d;
	    $best_d_error = $error;
	    $best = '*';
	};
	 print "  $best $trial/$j = $result  ($error)\n" if ($debug);
    };
    $last_a = $a;  $last_e = $e;
    $a = $c + $best_d;
    $i = $j;
};


exit 0;