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

Edited by mirod on 2003-02-13: added br/code/readmore tags
Edited by broquaint on 2003-02-14: turned URIs into links