Category: 
Miscellaneous 
Author/Contact Info 

Description: 
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
roundoff error.
Well the trick is called continued fractions
and the following program shows how you do it. In practice
you just do enough iterations that you get down to roundoff
error and stop with that.
I suggest trying out a few fractions, and after that some
expressions. Two fun expressions are
4 * atan2(1,1)
(1 + 5**0.5)/2
The first will give you some wellknown 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.
For more on the iterator techniques I used, I suggest
this
article.

#! /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;
};
}

Re: Continued Fractions by fundflow (Chaplain) on Nov 16, 2000 at 21:10 UTC 
Cool script. Thanks.
Small suggestion: It is usually useful to make scripts
pipeable, which is really easy for us perlers. All you need
to do is change the input to something like:
print "Please enter a number or expression: ";
my $num=<>  exit(0);
print "How many iterations: ";
my $count=<>  exit(0);
$num=eval($num);
chomp($count);
and then it will work from stdin and will exit on ^D (^Z on dos?)
if it reads from STDIN.
Cheers  [reply] [d/l] 

Why would <STDIN> differ from <> in this respect? Hitting an EOF character should work the same in either case. The only difference between <STDIN> and <> is that the latter will step through the lines of any files specified on the command line:
$ perl script.pl
$ perl script.pl file1 file2 filen
The first reads from the STDIN, but the second only references the contents of each of the files.  [reply] [d/l] [select] 

Yup, little difference indeed.
The emphasis was on the  exit thing.
Without it the script will have to be killed explicitly, and
perl script.pl < my_prepared_input > theoutput
will loop forever.
Anyway, it was just a small suggestion (and a good practice
in general)
 [reply] [d/l] [select] 

First of all as merlyn has pointed out, seeing interactive messages mixed with reading from @ARGV usually is a sign of problems.
Secondly I intentionally didn't make it pipeable. The point of this script is to let people figure out how it works, and then run it interactively (adding debug messages if you want) to give people a sense of how well it works. Any programmer should be able to figure out how to hit controlC. *shrug*
Before making it pipeable what I would do is add some logic to cut off at a sufficiently good approximation rather than trying to display the approximation process to the user. My not doing so is my way of saying that this code is not meant to be particularly useful in a pipeline. I didn't do that here because the entire point is education, not production. That is, rather than give a useful answer, it shows you the successive steps towards one.
In fact I probably would never make this pipeable. Instead it would go into a module, and you would use that module in your programs.
 [reply] 
Re: Continued Fractions by Dominus (Parson) on Nov 16, 2000 at 21:41 UTC 
 [reply] 

Simpler but far, far less reliable. I just tried it, and
tried some sample values. Specifically I entered (0.923076923076923) (ie 12/13) and it got into an infinite loop.
Didn't mean to do that. Meant to come back and point out that it didn't put this in an easily recognized form.
To make that point try 11.231 out. Mine will reduce that
to 11231/1000 in 5 steps while for me yours gets 160767181003805/14314591844342. (I am substantially less sensitive to roundoff errors until the iteration where I try to account for the roundoff error. Which is where you try to exit.)
 [reply] 

Hey, I only said it was simpler. I didn't say it was better.
:)
 [reply] 
(jeffa) Re: Continued Fractions by jeffa (Bishop) on Nov 16, 2000 at 22:54 UTC 
sub gcd {
my ($n, $m) = @_;
while ($m) {
my $k = $n % $m;
($n, $m) = ($m, $k);
}
return $n;
}
. . . particularly the use of the mod operator.
I remember reading Robert Sedgewick's
Mastering Algorithms in C++
and really getting a kick out of the first example  Euclide's
greatest common denominator formula.
At first, Sedgewick implemented it with the division
operator:
my $k = $n / $m;
Then he explained that replacing the / operator
with the % operator results in about 75% less iterations
through the while loop. Incredible.
Jeff  [reply] 
Back to Code Catacombs

