Symbolic algebra

[Brian Troutwine]

I'm trying to write a program to output a Lagrange interpolating
polynomial.
http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html

So far what I'm getting is output that looks somewhat along the lines of:

((x-2)(x-3)(x-4)(x-5)(1)/(24))+
((x-1)(x-3)(x-4)(x-5)(2)/(-6))+
((x-1)(x-2)(x-4)(x-5)(6)/(4))+
((x-1)(x-2)(x-3)(x-5)(24)/(-6))+
((x-1)(x-2)(x-3)(x-4)(120)/(24))

What I want is to be able to take that and output it in the form of

Ax^n + Bx^(n-1) ... + C

Does anybody have the experience to tell me how to calculate symbolically?
(Specifically how to implement a symbolic algebra module.)

 

[Michele Dondi]

I'm trying to write a program to output a Lagrange interpolating
polynomial.
[snip]

Does anybody have the experience to tell me how to calculate symbolically?
(Specifically how to implement a symbolic algebra module.)

There may well be some dedicated Perl module, but are you sure that
would be the best choice? There's a highly specialized and optimized
softare package called FORM especially aimed at tasks like yours (and
more complex, of course: it is particularly popular amongst high
energy physics researchers). Since it is compiler, after all, I think
you may efficiently use it from Perl if needed...

Once I did something similar myself, feeding it imput from Perl and
collecting its output with Perl to transform it into LaTeX code.
Unfortunately I don't have anything of that "work" left. It was really
nothing important however.


Michele

 

[M.J.T. Guy]

I'm trying to write a program to output a Lagrange interpolating
polynomial.
http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html

So far what I'm getting is output that looks somewhat along the lines of:

((x-2)(x-3)(x-4)(x-5)(1)/(24))+
((x-1)(x-3)(x-4)(x-5)(2)/(-6))+
((x-1)(x-2)(x-4)(x-5)(6)/(4))+
((x-1)(x-2)(x-3)(x-5)(24)/(-6))+
((x-1)(x-2)(x-3)(x-4)(120)/(24))

What I want is to be able to take that and output it in the form of

Ax^n + Bx^(n-1) ... + C

Implementing a general algebra package in Perl would be a substantial
task, but for your problem you need much less, so a full-blown algebra
package is massive overkill.

You only need to deal with polynomials in one variable (x), and only
need a few operations.

[ All following code untested ]

We represent a polynomial by a reference to an array of coefficients.

a) Create a constant polynomial $a
[ $a ];

b) Multiply the polynomial $b by (x-$a)
push @$b, 0; # degree increases by 1
for (my $i = $#$b; $i >= 0; $i--) { # NB must do backwards
$b->[$i] *= -$a;
$b->[$i] += $b->[$i-1] if $i; # don't reference $b[-1]
};

c) Add polynomial $a to $b
# next line not needed if your polynomials all have same degree
push @$b, (0) x (@a-@b) if @a > @b;
$b->[$_] += $a->[$_] for 0..$#a;

d) Stringify a polynomial $a
my $s = '';
for (my $i = $#$a; $i >= 0; $i--) {
my $c = $a->[$i] or next; # skip zero terms
$s .= ' ' if $s;
$s .= ($c > 0 ? "+ $c" : '- '.-$c);
$s .= 'x' if $i;
$s .= "^$i" if $i > 1;
};
print "The polynomial is $s\n";


Mike Guy

 

[Brian Troutwine]

Never did see a good answer to this one so I played around with it until
the following evolved. Uses the Math:olynomial module to simplify the
math - Gee! what a concept. Fred

use strict;
use warnings;
use diagnostics;
use Math:olynomial;

Math:olynomial->verbose(1);

my @Polys;
my $Total = Math:olynomial->new();
my @zeros = (-1,-2,-3,-4,-5);
my @weights = (1/24,-2/6,6/4,-24/6,120/24);

for my $p (0..4) {
my $Poly = Math:olynomial->new(0,1); for my $i (0..4) {
next if ($i == $p);
my $Term = Math:olynomial->new(1, $zeros[$i]); $Poly = $Poly *
$Term;
}
push @Polys, $Poly;
$Total += $Poly * $weights[$p];
}
$"="\n";
print "\nPolynomials are:\n@Polys\n\n$Total\n";

Neat method, but if you'll notice, the polynomial wasn't the same for
each term. Not a big problem as with a bit of tweaking yours could have
accounted for that (it's a very regular pattern, the trick is just
allowing for a potentially very large set of terms). Actually, I've
figured it out since last I checked here. And I managed to completely
avoid using computer algebra systems or even symbolic variables. Here's
most of my code (@zork is the set of numbers either piped in or read from
a file, as this is a program to analyze patterns from outside):


chomp @zork;

for( my $i = 0; $i <= $#zork; $i++){
my $flag = 0;
my $thisy = Math::BigInt->new($zork[$i]); my $denominator =
Math::BigInt->new(1); for( my $j = 0; $j <= $#zork; $j++){
if ($j != $i){
$denominator->bmul($i - $j);
}}
my @numerator = (0);
for( my $k = 0; $k <= $#zork; $k++){
if ($k != $i){
push @numerator,(-$k);
}}

print "$#numerator\n\n";

for(my $l =0; $l<=$#numerator; $l++){
if($l == 0){
@poly1 = (0,1);
$l++;
}

@poly2 = ($numerator[$l],1);

my @neweq = (0);

for(my $m=0; $m<=$#poly1; $m++){
for(my $n=0; $n<=$#poly2; $n++){
$neweq[$m+$n] += ($poly1[$m] * $poly2[$n]); ##############
}}
@poly1 = @neweq;
}
}
for(my $o=0; $o<=$#poly1; $o++){
my $coef = frac("$thisy","$denominator"); $poly1[$o] *= ($coef);
}
}
for(my $p=0; $p<=$#poly1; $p++){
$equation[$p] += $poly1[$p];
}
}

shift @equation;

my $flag = 0;

for(my $q = $#equation; $q >= 0; $q--){
unless($equation[$q] == 0){
if (($flag != 0) && ($equation[$q] > 0)){
print "+";
}
print $equation[$q];

unless($q == 0){
print "x";
}
if($q > 1){
print "^$q";
}
print "\n\n";
$flag++;
}}

if ($flag == 0){
print "0\n";
}


Essentially, instead of making a (x-5)... I set an array to (-5,...). Then
I make two poly arrays, where the zero position will have the x^0 value,
the one position the x^1, etc. This works very well, until we get up to
high numbers (around expanding 18 binomials).

(If you have any questions as to what certain parts do, please ask,
because it might mean I've got something useless that I haven't caught.)
Anyway, if you'll notice the line of #'s, that's where the program's
breaking at over 17 terms. If somebody could help me out here, I either
need to implement BigInt's or Fraction's and I'm having trouble.

Thanks.