Continued Fractions

[tilly]

Category:

Miscellaneous

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;
  };
}