A method for finding a polynomial of degree <= n-1 to fit n values given for the integers 1,2,3,...

Note that Math::Polynomial has its own version of this in the function "interpolate", which is probably faster. Also note that while polynomial interpolation has its uses, it's generally considered invalid to use it as I do here, since the interpolation is likely to be completely unrelated to the original function outside the bounds of the interpolation points.

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 w
+anted
# 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";
[download]

perl -pe'$\=($\*32768.5+ord)%4**8 .$/for/./gs}{'
[download]

What does this code do? It computes the same checksum as BSD's historic "sum" program. It's sum1 on http://ppt.perl.org/commands/sum/sum.theo