sourcecode
tilly
<code>
#! /usr/bin/perl -w
use strict;
while (1) {
print "Please enter a number or expression: ";
my $num = eval(scalar <STDIN>);
print "How many iterations: ";
chomp(my $count = <STDIN>);
print "Doing $count iterations of approximations to $num.\n";
my $f = ret_frac_iter($num);
for (1..$count) {
my ($n, $m) = $f->();
my $approx = $n/$m;
print "$n/$m = $approx\n";
}
print "\n";
}
# Takes a number, returns the best integer approximation and (in list
# context) the error.
sub best_int {
my $x = shift;
my $approx = sprintf '%.0f', $x;
if (wantarray) {
return ($approx, $x - $approx);
}
else {
return $approx;
}
}
# Takes a numerator and denominator, in scalar context returns
# the best fraction describing them, in list the numerator and
# denominator
sub frac_standard {
my $n = best_int(shift);
my $m = best_int(shift);
my $k = gcd($n, $m);
$n /= $k;
$m /= $k;
if ($m < 0) {
$n *= -1;
$m *= -1;
}
if (wantarray) {
return ($n, $m);
}
else {
return "$n/$m";
}
}
# Euclidean algorithm for calculating a GCD.
# Takes two integers, returns the greatest common divisor.
sub gcd {
my ($n, $m) = @_;
while ($m) {
my $k = $n % $m;
($n, $m) = ($m, $k);
}
return $n;
}
# Takes a list of terms in a continued fraction, and converts it
# into a fraction.
sub ints_to_frac {
my ($n, $m) = (0, 1); # Start with 0
while (@_) {
my $k = pop;
if ($n) {
# Want frac for $k + 1/($n/$m)
($n, $m) = frac_standard($k*$n + $m, $n);
}
else {
# Want $k
($n, $m) = frac_standard($k, 1);
}
}
return frac_standard($n, $m);
}
# Takes a number, returns an anon sub which iterates through a set of
# fractional approximations that converges very quickly to the number.
sub ret_frac_iter {
my $x = shift;
my $term_iter = ret_next_term_iter($x);
my @ints;
return sub {
push @ints, $term_iter->();
return ints_to_frac(@ints);
}
}
# terms of a continued fraction converging on that number.
sub ret_next_term_iter {
my $x = shift;
return sub {
(my $n, $x) = best_int($x);
if (0 != $x) {
$x = 1/$x;
}
return $n;
};
}
</code>
Everyone knows how to take a fraction and get a decimal
number out. What is not so obvious is how to take the
decimal and get a fraction back. Particularly with
round-off error.<P>
Well the trick is called <A HREF=http://home.flash.net/~mherk/contfrac.htm>continued fractions</A>
and the following program shows how you do it. In practice
you just do enough iterations that you get down to round-off
error and stop with that.<P>
I suggest trying out a few fractions, and after that some
expressions. Two fun expressions are
<code>
4 * atan2(1,1)
(1 + 5**0.5)/2
</code>
The first will give you some well-known approximations of
pi, and the second is the golden mean. The golden mean has
the distinction of being the slowest continued fraction to
converge. Its terms are also interesting to look at.<P>
For more on the iterator techniques I used, I suggest
<A HREF=http://perl.plover.com/Stream/stream.html>this
article</A>.
Miscellaneous