Symbolic algebra

Discussion in 'Perl Misc' started by Brian Troutwine, Sep 30, 2004.

  1. 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.)
    Brian Troutwine, Sep 30, 2004
    #1
    1. Advertising

  2. On Thu, 30 Sep 2004 14:59:35 -0500, Brian Troutwine
    <> wrote:

    >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
    --
    {$_=pack'B8'x25,unpack'A8'x32,$a^=sub{pop^pop}->(map substr
    (($a||=join'',map--$|x$_,(unpack'w',unpack'u','G^<R<Y]*YB='
    ..'KYU;*EVH[.FHF2W+#"\Z*5TI/ER<Z`S(G.DZZ9OX0Z')=~/./g)x2,$_,
    256),7,249);s/[^\w,]/ /g;$ \=/^J/?$/:"\r";print,redo}#JAPH,
    Michele Dondi, Oct 1, 2004
    #2
    1. Advertising

  3. Brian Troutwine

    M.J.T. Guy Guest

    Brian Troutwine <> wrote:
    >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
    M.J.T. Guy, Oct 5, 2004
    #3
  4. On Tue, 05 Oct 2004 23:41:34 +0000, Fred Toewe wrote:

    > Never did see a good answer to this one so I played around with it until
    > the following evolved. Uses the Math::polynomial module to simplify the
    > math - Gee! what a concept. Fred
    >
    > use strict;
    > use warnings;
    > use diagnostics;
    > use Math::polynomial;
    >
    > Math::polynomial->verbose(1);
    >
    > my @Polys;
    > my $Total = Math::polynomial->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::polynomial->new(0,1); for my $i (0..4) {
    > next if ($i == $p);
    > my $Term = Math::polynomial->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.
    Brian Troutwine, Oct 7, 2004
    #4
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. ckumar
    Replies:
    2
    Views:
    447
    ckumar
    Jan 17, 2005
  2. Robert M. Gary

    Matrix Algebra support

    Robert M. Gary, Nov 11, 2005, in forum: Java
    Replies:
    4
    Views:
    451
    Googmeister
    Nov 12, 2005
  3. Replies:
    6
    Views:
    5,076
  4. xhunga

    symbolic algebra computation with list

    xhunga, Apr 27, 2007, in forum: C Programming
    Replies:
    1
    Views:
    306
  5. Brian Troutwine

    Symbolic algebra

    Brian Troutwine, Sep 30, 2004, in forum: Perl Misc
    Replies:
    0
    Views:
    69
    Brian Troutwine
    Sep 30, 2004
Loading...

Share This Page