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