http://qs321.pair.com?node_id=41961
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 round-off 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 round-off 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 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.

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;
  };
}
Replies are listed 'Best First'.
(jeffa) Re: Continued Fractions
by jeffa (Bishop) on Nov 16, 2000 at 22:54 UTC
    You rule Tilly!

    I had to point out the mod operator in this function:

    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

Re: Continued Fractions
by Dominus (Parson) on Nov 16, 2000 at 21:41 UTC
      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.)

        Hey, I only said it was simpler. I didn't say it was better.
        :)
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

      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 control-C. *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.

      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.
        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)