fizbin's scratchpad

Public Scratchpad

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]

Regular expression speed comparison, or a meditation on such
I still need to add bleadperl numbers to this, and I may be redoing all my numbers soon enough (I'm getting a spiffier machine), but 5.9.3 was added and, well, it didn't do as well as it should have.

Note to self: Remember that Tilly came up with this in a chatterbox golfing session:

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