darek56

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

Posted to the free world!!!
No warranties nor guaranties, free to use and enjoy.
Bidek

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

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

use Data:umper;

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

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