# Lagrange Interpolating Polynomial

#### Do you have a question? Post it now! No Registration Necessary.  Now with pictures!

•  Subject
• Author
• Posted on
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.

## Re: Lagrange Interpolating Polynomial

Brian Troutwine wrote:

[ snip ]

\$_ *= \$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/">http://johnbokma.com/mexit/
personal page:       http://johnbokma.com/">http://johnbokma.com/
Experienced programmer available:     http://castleamber.com/">http://castleamber.com/
Happy Customers: http://castleamber.com/testimonials.html">http://castleamber.com/testimonials.html

## Re: Lagrange Interpolating Polynomial

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

No warnings?

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

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

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

The for-loop can be written:

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

Similarly:

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

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

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;

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

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

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

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.

# 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

## Re: Lagrange Interpolating Polynomial

Anno Siegel wrote:

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

## Re: Lagrange Interpolating Polynomial

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

Anno

## Re: Lagrange Interpolating Polynomial

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.
ben@morrow.me.uk                                    The Levellers, 'Believers'

## Re: Lagrange Interpolating Polynomial

On Tue, 12 Oct 2004 12:48:52 +0100, Ben Morrow wrote:

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.)

## Re: Lagrange Interpolating Polynomial

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,                * ben@morrow.me.uk
it is worth reading when you are sixty." - C.S.Lewis

## Re: Lagrange Interpolating Polynomial

On Wed, 13 Oct 2004 21:13:53 +0100, Ben Morrow wrote:

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.

## Re: Lagrange Interpolating Polynomial

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.me.uk

## Re: Lagrange Interpolating Polynomial

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

## Re: Lagrange Interpolating Polynomial

[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

## Re: Lagrange Interpolating Polynomial

On 13 Oct 2004 07:36:19 GMT, anno4000@lublin.zrz.tu-berlin.de (Anno
Siegel) wrote:

<OT>

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

No, no, it was programmed in Forth.  See Genesis 1:12:
"And the earth brought Forth ..."
- Robert Israel in sci.math, thread "Why numbers?"

No, I'm not.

Anno

## Re: Lagrange Interpolating Polynomial

On 13 Oct 2004 16:26:30 GMT, anno4000@lublin.zrz.tu-berlin.de (Anno
Siegel) wrote:

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

Michele
--
->(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,

## Re: Lagrange Interpolating Polynomial

On Tue, 12 Oct 2004 19:48:29 -0700, Bidek wrote:

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.

## Re: Lagrange Interpolating Polynomial

[A complimentary Cc of this posting was sent to
Brian Troutwine

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

Hope this helps,
Ilya