Lagrange Interpolating Polynomial

Discussion in 'Perl Misc' started by Brian Troutwine, Oct 12, 2004.

  1. Here's my finished Lagrange polynomial finder. It runs, and unless a
    stress test suggests otherwise (to be performed overnight), works
    excellently, but merely is really slow. Any suggestions to get it running
    faster?

    #!/usr/bin/perl
    use diagnostics;
    use strict;
    use Getopt::Std;
    use Math::BigRat;

    #Variables
    #zork is a common array to allow for different input sources
    #The starting zero allows us to use this algorithm for as few
    #as two input values. Hopefully the next version will allow
    #for one term, but that's trivial to do by hand.
    my @zork = (0);
    our @equation = (0);
    our @poly1 = ("");
    our @poly2 = ("");
    our $opt_f = "";
    getopt("f");

    #Opens pipe if no -f used
    if($opt_f eq ""){
    foreach my $foo (<>){
    push @zork, Math::BigRat->new($foo);
    }}

    #Opens file specified by -f
    else{
    open(FOO, "<$opt_f") || die("File does not exist.");
    my @bar=<FOO>;
    close(FOO);
    foreach my $bar (@bar){
    push @zork, Math::BigRat->new($bar);
    }}

    chomp @zork;

    #$i represents the number of the term we're working with.
    for( my $i = 0; $i <= $#zork; $i++){
    my $thisy = Math::BigRat->new($zork[$i]);

    #Creates denominator for each term.
    my $denominator = Math::BigRat->new('1/1');
    for( my $j = 0; $j <= $#zork; $j++){
    if ($j != $i){
    $denominator *= ($i - $j);
    }}

    #Makes array of n's in (x-n) where n is
    #all integers from 0 to $#zork except
    #$i.
    my @numerator;
    for( my $k = 0; $k <= $#zork; $k++){
    if ($k != $i){
    push @numerator,(-$k);
    }}

    #Computes numerator as polynomial
    for(my $l = 0; $l<=$#numerator; $l++){

    #The polynomials will be arrays in which any position
    #will hold the coefficient of the corresponding power
    #of x. For example, the first position (0) will hold
    #the constant (the coefficient for x^0).
    if($l == 0){
    @poly1 = (0,1);
    $l++;
    }

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

    my @neweq = (0);

    #To multiply our polynomials, we take each of the
    #terms of the first polynomial, and multiply them
    #by each of the terms of the second, and add the
    #resulting coefficient to the sum of their positions
    for(my $m=0; $m<=$#poly1; $m++){
    for(my $n=0; $n<=$#poly2; $n++){
    $neweq[$m+$n] += Math::BigRat->new($poly1[$m] * $poly2[$n]);
    }}

    #Set poly1 to the new equation for the next loop through
    @poly1 = @neweq;
    }

    #To avoid decimals, we've used fractions to multiply
    #each of our terms by.
    for(my $o=0; $o<=$#poly1; $o++){
    $poly1[$o] *= $thisy / $denominator;
    }

    #Adding up our polynomials, because we want to simplify
    for(my $p=0; $p<=$#poly1; $p++){
    $equation[$p] += $poly1[$p];
    }
    }

    my $flag = 0;

    #Output
    for(my $q = $#equation; $q >= 0; $q--){

    #Doesn't output zero terms
    unless($equation[$q] == 0){

    #No plus if it's first or negative, it just makes sense
    if (($flag != 0) && ($equation[$q] > 0)){
    print "+";
    }
    print $equation[$q];

    #Prints x except for on constant, and power except for
    #first power
    unless($q == 0){
    print "x";
    }
    if($q > 1){
    print "^$q";
    }
    print "\n\n";
    $flag++;
    }}

    #In case somebody happens to only input zeroes, we don't want
    #to leave that poor bastard with no output.
    if ($flag == 0){
    print "0\n";
    }


    Thanks.
     
    Brian Troutwine, Oct 12, 2004
    #1
    1. Advertising

  2. Brian Troutwine

    John Bokma Guest

    Brian Troutwine wrote:

    > Here's my finished Lagrange polynomial finder. It runs, and unless a
    > stress test suggests otherwise (to be performed overnight), works
    > excellently, but merely is really slow. Any suggestions to get it running
    > faster?


    [ snip ]

    > for(my $o=0; $o<=$#poly1; $o++){
    > $poly1[$o] *= $thisy / $denominator;
    > }


    $_ *= $thisy / $denominator for @poly1;

    ?

    I see a lot of indexed looping, probably can be rewritten as the one above.


    --
    John MexIT: http://johnbokma.com/mexit/
    personal page: http://johnbokma.com/
    Experienced programmer available: http://castleamber.com/
    Happy Customers: http://castleamber.com/testimonials.html
     
    John Bokma, Oct 12, 2004
    #2
    1. Advertising

  3. Brian Troutwine

    Anno Siegel Guest

    Brian Troutwine <> wrote in comp.lang.perl.misc:
    > Here's my finished Lagrange polynomial finder. It runs, and unless a
    > stress test suggests otherwise (to be performed overnight), works
    > excellently, but merely is really slow. Any suggestions to get it running
    > faster?


    Not significantly. You have a few unnecessary calls to
    Math::BigRat->new, but otherwise you seem to do what must be done.

    > #!/usr/bin/perl
    > use diagnostics;
    > use strict;
    > use Getopt::Std;
    > use Math::BigRat;


    No warnings?

    > #Variables
    > #zork is a common array to allow for different input sources
    > #The starting zero allows us to use this algorithm for as few
    > #as two input values. Hopefully the next version will allow
    > #for one term, but that's trivial to do by hand.
    > my @zork = (0);
    > our @equation = (0);
    > our @poly1 = ("");
    > our @poly2 = ("");
    > our $opt_f = "";
    > getopt("f");


    It is preferable to declare variables immediately before their first
    use.

    > #Opens pipe if no -f used
    > if($opt_f eq ""){
    > foreach my $foo (<>){
    > push @zork, Math::BigRat->new($foo);
    > }}
    >
    > #Opens file specified by -f
    > else{
    > open(FOO, "<$opt_f") || die("File does not exist.");
    > my @bar=<FOO>;
    > close(FOO);
    > foreach my $bar (@bar){
    > push @zork, Math::BigRat->new($bar);
    > }}
    >
    > chomp @zork;


    @zork is an array of BigRat's. Chomping makes no sense.

    > #$i represents the number of the term we're working with.
    > for( my $i = 0; $i <= $#zork; $i++){
    > my $thisy = Math::BigRat->new($zork[$i]);


    $zork[ $i] is already a BigRat object. No need for ->new.

    > #Creates denominator for each term.
    > my $denominator = Math::BigRat->new('1/1');
    > for( my $j = 0; $j <= $#zork; $j++){
    > if ($j != $i){
    > $denominator *= ($i - $j);
    > }}


    The for-loop can be written:

    $denominator *= ( $i - $_) for grep $_ != $i, 0 .. $#zork;

    >
    > #Makes array of n's in (x-n) where n is
    > #all integers from 0 to $#zork except
    > #$i.
    > my @numerator;
    > for( my $k = 0; $k <= $#zork; $k++){
    > if ($k != $i){
    > push @numerator,(-$k);
    > }}


    Similarly:

    my @numerator = map -$_, grep $_ != $i, 0 .. $#zork;

    > #Computes numerator as polynomial
    > for(my $l = 0; $l<=$#numerator; $l++){
    >
    > #The polynomials will be arrays in which any position
    > #will hold the coefficient of the corresponding power
    > #of x. For example, the first position (0) will hold
    > #the constant (the coefficient for x^0).
    > if($l == 0){
    > @poly1 = (0,1);
    > $l++;
    > }


    Doing this inside the loop is awkward. Initialize $poly1 before
    the loop, and start at $l = 1.

    > @poly2 = ($numerator[$l],1);
    >
    > my @neweq = (0);
    >
    > #To multiply our polynomials, we take each of the
    > #terms of the first polynomial, and multiply them
    > #by each of the terms of the second, and add the
    > #resulting coefficient to the sum of their positions
    > for(my $m=0; $m<=$#poly1; $m++){
    > for(my $n=0; $n<=$#poly2; $n++){
    > $neweq[$m+$n] += Math::BigRat->new($poly1[$m] * $poly2[$n]);
    > }}


    You are using a full-fledged multiplication routine (which, in a
    greater scheme of things, would deserve to be a subroutine). But
    all you do is multiply a polynomial with a linear expression (x - a).
    That could be done in an ad-hoc manner thusly:

    my @neweq = map $_ * $numerator[ $l], @poly1;
    $neweq[ $_ + 1] += $poly1[ $_] for 0 .. $#poly1;

    > #Set poly1 to the new equation for the next loop through
    > @poly1 = @neweq;
    > }
    >
    > #To avoid decimals, we've used fractions to multiply
    > #each of our terms by.
    > for(my $o=0; $o<=$#poly1; $o++){
    > $poly1[$o] *= $thisy / $denominator;
    > }


    $_ *= $thisy / $denominator for @poly1;

    > #Adding up our polynomials, because we want to simplify
    > for(my $p=0; $p<=$#poly1; $p++){
    > $equation[$p] += $poly1[$p];
    > }


    $equation[ $_] += $poly1[ $_] for 0 .. $#poly1;

    > }
    >
    > my $flag = 0;


    That belongs to the output routine and should be below the
    comment.

    > #Output


    I'd use an auxiliary routine to format one term for the output.
    Quite generally, your code could gain a lot if you used more
    subroutines to factor out common operations.

    > for(my $q = $#equation; $q >= 0; $q--){
    >
    > #Doesn't output zero terms
    > unless($equation[$q] == 0){
    >
    > #No plus if it's first or negative, it just makes sense
    > if (($flag != 0) && ($equation[$q] > 0)){
    > print "+";
    > }
    > print $equation[$q];
    >
    > #Prints x except for on constant, and power except for
    > #first power
    > unless($q == 0){
    > print "x";
    > }
    > if($q > 1){
    > print "^$q";
    > }
    > print "\n\n";
    > $flag++;
    > }}
    >
    > #In case somebody happens to only input zeroes, we don't want
    > #to leave that poor bastard with no output.
    > if ($flag == 0){
    > print "0\n";
    > }


    # collect non-zero terms
    my @terms = grep $_->[ 0], map [ $equation[ $_], $_], 0 .. $#equation;

    my $str = '';
    $str .= term_to_str( @$_) for reverse @terms;
    $str =~ s/^\+//;
    print $str || '0';

    sub term_to_str {
    my ( $fac, $exp) = @_;
    my $power = '';
    $power .= 'x' if $exp;
    $power .= "^$exp" if $exp != 1;
    $fac = "+$fac" unless $fac < 0;
    "$fac$power\n\n";
    }

    Anno
     
    Anno Siegel, Oct 12, 2004
    #3
  4. Brian Troutwine

    Ben Morrow Guest

    Quoth Brian Troutwine <>:
    > Here's my finished Lagrange polynomial finder. It runs, and unless a
    > stress test suggests otherwise (to be performed overnight), works
    > excellently, but merely is really slow. Any suggestions to get it running
    > faster?
    >
    > #!/usr/bin/perl
    > use diagnostics;
    > use strict;
    > use Getopt::Std;
    > use Math::BigRat;


    I'm not competent to comment on the algorithm (and anyway Anno already
    has) but I would suspect the default BigRat implementation is the
    slowest. You could try changing to ::GMP or something to see if it makes
    a difference.

    Ben

    --
    If you put all the prophets, | You'd have so much more reason
    Mystics and saints | Than ever was born
    In one room together, | Out of all of the conflicts of time.
    The Levellers, 'Believers'
     
    Ben Morrow, Oct 12, 2004
    #4
  5. Anno Siegel wrote:

    > Brian Troutwine <> wrote in comp.lang.perl.misc:
    >> Here's my finished Lagrange polynomial finder. It runs, and unless a
    >> stress test suggests otherwise (to be performed overnight), works
    >> excellently, but merely is really slow. Any suggestions to get it running
    >> faster?

    >
    > Not significantly. You have a few unnecessary calls to
    > Math::BigRat->new, but otherwise you seem to do what must be done.
    >
    >> #!/usr/bin/perl
    >> use diagnostics;
    >> use strict;
    >> use Getopt::Std;
    >> use Math::BigRat;

    >
    > No warnings?
    >

    He's got use diagnostics. Use warnings would be redundant.

    --
    Christopher Mattern

    "Which one you figure tracked us?"
    "The ugly one, sir."
    "...Could you be more specific?"
     
    Chris Mattern, Oct 12, 2004
    #5
  6. Brian Troutwine

    Anno Siegel Guest

    Chris Mattern <> wrote in comp.lang.perl.misc:
    > Anno Siegel wrote:
    >
    > > Brian Troutwine <> wrote in comp.lang.perl.misc:
    > >> #!/usr/bin/perl
    > >> use diagnostics;
    > >> use strict;
    > >> use Getopt::Std;
    > >> use Math::BigRat;

    > >
    > > No warnings?
    > >

    > He's got use diagnostics. Use warnings would be redundant.


    Oh, right. Well, that adds some sluggishness too.

    Anno
     
    Anno Siegel, Oct 12, 2004
    #6
  7. On Tue, 12 Oct 2004 12:48:52 +0100, Ben Morrow wrote:

    > You could try changing to ::GMP or something to see if it makes
    > a difference.


    Unfortunately, I need rational number support, and since Math::Fraction
    doesn't work particularly well, BigRat seems my best option. Do you know
    of any better? (I looked into ::GMP, but that seemed to be large numbers
    and not large rationals numbers.)
     
    Brian Troutwine, Oct 13, 2004
    #7
  8. Brian Troutwine

    Bidek Guest

    Brian,
    Have you tries the Newton's method instead, here's a suggestion.

    The following Perl program estimates the roots of a polynomial using
    Newton's method and Horner's rule.

    #################################################

    #!/opt/perl/5.8.0/bin/perl -w

    use Data::Dumper;

    use constant EPSILON => 10 ** -10;
    use constant { DEBUG => 0, OUTPUT => 1 };

    my @polynomial = qw( 4 -44 165 -224 49 );

    print_poly (@polynomial);
    @roots = NewtonMethod (@polynomial);
    print_roots (@roots);

    exit 0;

    sub NewtonMethod {

    local (@P) = @_;
    local (@Q, @PP, @QP, @results) = @P;
    local ($count) = ($#P);

    # calculate the derivative of the polynomial
    for ($i=1; $i <= $count; $i++) {
    $QP[$i - 1] = $PP[$i - 1] = $i * $P[$i];
    }

    do {
    # initialize the root with a "small" random number
    $x = rand(); # INIT_X;
    $val = HornerRule ($count - 1, $x, @Q);

    $iteration = 0;

    print "=" x 40, "\n" if OUTPUT;

    LOOP:
    $iteration ++;

    my $dx = $val / HornerRule ($count - 1, $x, @QP);

    if (abs($dx) > EPSILON) {
    for ($c = 0; $c < 40; $c++) {
    $y = $x - $dx;
    $tmp_val = HornerRule($count, $y, @Q);

    printf "I: %3d X: %13.10f Px: %13.10f\n", $iteration, $y, $tmp_val if OUTPUT;

    if ($tmp_val ** 2 < $val ** 2) {
    $x = $y, $val = $tmp_val;
    goto LOOP;
    }
    $dx /= 4.0;
    }
    }

    # evaluate the polynomial and its derivative
    $Px = HornerRule($#P, $x, @P);
    $PPx = HornerRule($#P - 1, $x, @PP);

    # apply Newton's formula
    if ($PPx != 0.0) {
    $x -= $Px / $PPx;
    $Px = HornerRule($#P, $x, @P);
    $PPx = HornerRule($#P - 1, $x, @PP);

    # apply Newton's formula a second time
    $x -= $Px / $PPx if $PPx;
    }

    if (abs($x) > EPSILON) {

    push @results, $x;

    # deflate the polynomial by a linear factor
    $A[$count - 1] = $Q[$count];
    for ($i = $count - 1; $i >= 1; $i--) {
    $A[$i - 1] = $Q[$i] + $x * $A[$i];
    }
    $count--;
    @Q = @A;
    }

    # calculate the derivative of Q(x)
    for ($i = 1; $i <= $count; $i++) {
    $QP[$i - 1] = $i * $Q[$i];
    }

    } while $count > 0;

    print Dumper @results if DEBUG;

    @results;
    }

    sub HornerRule {

    local ($n, $x, @P) = @_;

    local ($result) = $P[$n];

    print "=" x 20,"\n", Dumper $n, $x, @P if DEBUG;

    for ($i=$n-1; $i >= 0; $i--) {
    $result = $result * $x + $P[$i];
    }

    $result;
    }

    sub print_roots {

    local (@R) = @_;

    print "\nThe roots are: \n" , '-' x 20, "\n";
    for ($i=0; $i<=$#R; $i++) {
    printf "%d = %16.13f\n", $i+1, $R[$i];
    }
    print "-" x 20,"\n";
    }

    sub print_poly {

    local (@P) = @_;

    print "-" x 20,"\n";
    print "The polynomial is: \n", '-' x 20, "\n";
    for ($i = $#P; $i>0; $i--) {
    printf "%5d + 0.0x * X^%d + \n", $P[$i], $i;
    }

    printf "%5d + 0.0x\n", $P[0], $i;
    print "-" x 20,"\n";

    return 0;
    }
     
    Bidek, Oct 13, 2004
    #8
  9. Brian Troutwine

    Anno Siegel Guest

    Bidek <> wrote in comp.lang.perl.misc:
    > Brian,
    > Have you tries the Newton's method instead, here's a suggestion.
    >
    > The following Perl program estimates the roots of a polynomial using
    > Newton's method and Horner's rule.


    [code snipped]

    That's zero-finding, not interpolation -- an entirely different
    problem.

    There *is* a Newtonian method (different from Lagrange's) for polynomial
    interpolation, but that doesn't seem to be what you are talking about.

    Anno
     
    Anno Siegel, Oct 13, 2004
    #9
  10. On 13 Oct 2004 07:36:19 GMT, -berlin.de (Anno
    Siegel) wrote:

    <OT>
    >There *is* a Newtonian method (different from Lagrange's) for polynomial
    >interpolation, but that doesn't seem to be what you are talking about.


    Since we're taliking about mathematics now, by any chance do you
    happen to be a relative to C.L. Siegel?

    (I would have sent you an e-mail, but I don't know if you accept
    unsolicited correspondence.)
    </OT>


    Michele
    --
    >It's because the universe was programmed in C++.

    No, no, it was programmed in Forth. See Genesis 1:12:
    "And the earth brought Forth ..."
    - Robert Israel in sci.math, thread "Why numbers?"
     
    Michele Dondi, Oct 13, 2004
    #10
  11. On Tue, 12 Oct 2004 19:48:29 -0700, Bidek wrote:

    > Brian,
    > Have you tries the Newton's method instead, here's a suggestion.


    Thanks, but I'm not looking for roots, I'm fitting a curve to a set of
    points. I've got my algorithm down, I'm just trying to speed things up a
    bit (using BigRat really kills the time).

    On the other hand, I'll probably end up writing more pattern finding
    programs, and if I need roots, I'll give yours a shot.

    Thanks again.
     
    Brian Troutwine, Oct 13, 2004
    #11
  12. Brian Troutwine

    Anno Siegel Guest

    Michele Dondi <> wrote in comp.lang.perl.misc:
    > On 13 Oct 2004 07:36:19 GMT, -berlin.de (Anno
    > Siegel) wrote:
    >
    > <OT>
    > >There *is* a Newtonian method (different from Lagrange's) for polynomial
    > >interpolation, but that doesn't seem to be what you are talking about.

    >
    > Since we're taliking about mathematics now, by any chance do you
    > happen to be a relative to C.L. Siegel?


    No, I'm not.

    Anno
     
    Anno Siegel, Oct 13, 2004
    #12
  13. Brian Troutwine

    Ben Morrow Guest

    Quoth Brian Troutwine <>:
    > On Tue, 12 Oct 2004 12:48:52 +0100, Ben Morrow wrote:
    >
    > > You could try changing to ::GMP or something to see if it makes
    > > a difference.

    >
    > Unfortunately, I need rational number support, and since Math::Fraction
    > doesn't work particularly well, BigRat seems my best option. Do you know
    > of any better? (I looked into ::GMP, but that seemed to be large numbers
    > and not large rationals numbers.)


    I know only what it says in the Math::BigRat (0.10) perldoc, which
    implies that Math::BigRat uses Math::BigInt ('s internals) to do all its
    calculations, and that

    use Math::BigRat lib => 'GMP';

    will use Math::BigInt::GMP to do your BigRat stuff (if it's installed).

    Ben

    --
    "If a book is worth reading when you are six, *
    it is worth reading when you are sixty." - C.S.Lewis
     
    Ben Morrow, Oct 13, 2004
    #13
  14. On Wed, 13 Oct 2004 21:13:53 +0100, Ben Morrow wrote:
    > use Math::BigRat lib => 'GMP';
    >
    > will use Math::BigInt::GMP to do your BigRat stuff (if it's installed).
    >
    > Ben


    Status report:

    Well, It would appear that I have it since diagnostics doesn't spit me out
    any warning, but any speed advantage it has is negligible, I've tested
    multiple times, sometimes with yours coming out faster, sometimes with the
    original.

    Thanks anyway.
     
    Brian Troutwine, Oct 14, 2004
    #14
  15. Brian Troutwine

    Ben Morrow Guest

    Quoth Brian Troutwine <>:
    > On Wed, 13 Oct 2004 21:13:53 +0100, Ben Morrow wrote:
    > > use Math::BigRat lib => 'GMP';
    > >
    > > will use Math::BigInt::GMP to do your BigRat stuff (if it's installed).
    > >
    > > Ben

    >
    > Status report:
    >
    > Well, It would appear that I have it since diagnostics doesn't spit me out
    > any warning, but any speed advantage it has is negligible, I've tested
    > multiple times, sometimes with yours coming out faster, sometimes with the
    > original.


    That's (maybe) because you don't have it :)

    Try perl -MMath::BigInt::GMP -e1. If that succeeds, then you *do* have
    it. Math::BigRat will automatically and silently fall back to the
    pure-perl implementation if the one you ask for doesn't exist.

    Ben

    --
    We do not stop playing because we grow old;
    we grow old because we stop playing.
     
    Ben Morrow, Oct 14, 2004
    #15
  16. On 13 Oct 2004 16:26:30 GMT, -berlin.de (Anno
    Siegel) wrote:

    >> Since we're taliking about mathematics now, by any chance do you
    >> happen to be a relative to C.L. Siegel?

    >
    >No, I'm not.


    Well, you could say your cds are Siegel discs, anyway! ;-)


    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 14, 2004
    #16
  17. [A complimentary Cc of this posting was sent to
    Brian Troutwine
    <>], who wrote in article <>:
    > use Math::BigRat;


    Did you try Math::pari? It may even have interpolation builtin...
    Right, it is polint()...

    Hope this helps,
    Ilya
     
    Ilya Zakharevich, Oct 14, 2004
    #17
    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. Chen L.
    Replies:
    2
    Views:
    912
    Dik T. Winter
    Jul 6, 2004
  2. jitasi
    Replies:
    1
    Views:
    779
    Terry Reedy
    Mar 4, 2007
  3. Replies:
    0
    Views:
    842
  4. John Elrick
    Replies:
    9
    Views:
    122
    John Carter
    Apr 7, 2005
  5. Lloyd Zusman
    Replies:
    4
    Views:
    222
    Ken Bloom
    Mar 18, 2009
Loading...

Share This Page