use bigrat;
use strict;

# These are library routines for polynomial arithmetic - I could have used
# Math::Polynomial, but it's incompatible with "use bigrat", which I wanted
# to use.

# A polynomial is represented as an arrayref - the first term is the
# constant term, the next term is the x^1 term, then x^2, etc.

sub polymult {
  my ($a, $b) = @_;
  if (ref($a) ne 'ARRAY') {$a = [ $a ];}
  if (ref($b) ne 'ARRAY') {$b = [ $b ];}
  my $ans = [];
  for my $i (0..$#{$b}) {
    $ans = polyadd($ans, [ (0) x $i, map {$_*($b->[$i])} @$a ]);
  }
  $ans;
}

sub polyadd {
  my ($a, $b) = @_;
  if (ref($a) ne 'ARRAY') {$a = [ $a ];}
  if (ref($b) ne 'ARRAY') {$b = [ $b ];}
  if ($#{$a} < $#{$b}) {
    ($a, $b) = ($b, $a);
  }
  my $ans = [ @$a ];
  for my $i (0..$#{$b}) {
    $ans->[$i] += $b->[$i];
  }
  $ans;
}

sub polyeval {
  my ($x, $a) = @_;
  
  if (ref($a) != 'ARRAY') {return $a;}
  if (@$a == 0) {return 0;}
  if (@$a == 1) {return $a->[0];}
  my $prev = polyeval($x, [ @{$a}[1..$#$a] ]);
  my $ans = ($a->[0] + $x * $prev);
  $ans;
}

# Now on with the show

# Here place the given numbers plus whatever other ones you want
# to have show up later
my @values = qw ( 1 2 3 4 29 52 125 );

# The values you want to fit
my @xvals = (1..scalar(@values));

# The values you want to check
my @exvals = (@xvals, scalar(@values)+1);

my $P = [ shift @values ];

my $i = 0;

my $Q = 1;

while (@values) {
  $Q = polymult($Q , [-$xvals[$i],1]);
  $i++;
  my $newval = shift @values;
  my $oldval = polyeval($xvals[$i], $P);
  my $ratio  = polyeval($xvals[$i], $Q);
  
  next if ($newval == $oldval);
  my $addition = polymult($Q, ($newval - $oldval) / $ratio);
  $P = polyadd($P, $addition);
}

print "Polynomial: [ @$P ]\n";
print "cross-check of series, plus one more term:\n";
for my $v (@exvals) {print polyeval($v, $P), " ";}
print "\n";