/ 
Math-Polynomial-Solve-2.60
River stage zero No dependents
/ Math::Polynomial::Solve

NAME

Math::Polynomial::Solve - Find the roots of polynomial equations.

SYNOPSIS

use Math::Complex;  # The roots may be complex numbers.
use Math::Polynomial::Solve qw(poly_roots);

my @x = poly_roots(@coefficients);

or

use Math::Complex;  # The roots may be complex numbers.
use Math::Polynomial::Solve qw(:numeric);  # See the EXPORT section

#
# Find roots using the matrix method.
#
my @x = poly_roots(@coefficients);

#
# Alternatively, use the classical methods instead of the matrix
# method if the polynomial degree is less than five.
#
poly_option(hessenberg => 0);
@x = poly_roots(@coefficients);

or

use Math::Complex;  # The roots may be complex numbers.
use Math::Polynomial::Solve qw(:classical);  # See the EXPORT section

#
# Find the polynomial roots using the classical methods.
#

# Find the roots of ax + b
my @x1 = linear_roots($a, $b);

# Find the roots of ax**2 + bx +c
my @x2 = quadratic_roots($a, $b, $c);

# Find the roots of ax**3 + bx**2 +cx + d
my @x3 = cubic_roots($a, $b, $c, $d);

# Find the roots of ax**4 + bx**3 +cx**2 + dx + e
my @x4 = quartic_roots($a, $b, $c, $d, $e);

or

use Math::Polynomial::Solve qw(:utility);

my @coefficients = (89, 23, 23, 432, 27);

# Return a version of the polynomial with no leading zeroes
# and the leading coefficient equal to 1.
my @monic_form = simplified_form(@coefficients);

# Find the y-values of the polynomial at selected x-values.
my @xvals = (0, 1, 2, 3, 5, 7);
my @yvals = poly_evaluate(\@coefficients, \@xvals);

or

use Math::Polynomial::Solve qw(:sturm);

# Find the number of unique real roots of the polynomial.
my $no_of_unique_roots = poly_real_root_count(@coefficients);

DESCRIPTION

This package supplies a set of functions that find the roots of polynomials, along with some utility functions.

Functions making use of the Sturm sequence are also available, letting you find the number of real roots present in a range of X values.

In addition to the root-finding functions, the internal functions have also been exported for your use.

Numeric Functions

These are the functions that calculate the polynomial's roots through numeric algorithms. They are all exported under the tag "numeric".

poly_option()

Set options that affect the behavior of the poly_roots() function. All options are set to either 1 ("on") or 0 ("off").

This is the option function that deprecates set_hessenberg() and get_hessenberg().

Options may be set and saved:

#
# Set a few options...
#
poly_option(hessenberg => 0, root_function => 1);

#
# Get all of the current options and their values.
#
my %all_options = poly_option();

#
# Set some options but save the former option values
# for later.
#
my %changed_options = poly_option(hessenberg => 1, varsubst => 1);

The current options available are:

hessenberg

Use the QR Hessenberg matrix method to solve the polynomial. By default, this is set to 1. If set to 0, poly_roots() uses one of the classical root-finding functions listed below, if the degree of the polynomial is four or less.

root_function

If the polynomial is of the form ax**n + c, poly_roots() will make use of the root() function from Math::Complex. This takes precedence over the other solving methods.

varsubst

Reduce polynomials to a lower degree through variable substitution, if possible.

For example, with varsubst set to one and the polynomial to solve is 9x**6 + 128x**3 + 21, then poly_roots() will reduce the polynomial to 9y**2 + 128y + 21, where y = x**3, and solve the quadratic (either classically or numerically, depending on the hessenberg option). Taking the cube root of each quadratic root completes the operation.

This has the benefit of having a simpler matrix to solve, or if the hessenberg option is set to zero, has the effect of being able to use one of the classical methods on a polynomial of high degree. In the above example, the order-six polynomial gets solved with the quadratic_roots() function if the hessenberg option is zero.

Currently the variable substitution is fairly simple and will only look for gaps of zeros in the coefficients that are multiples of the prime numbers less than or equal to 31 (2, 3, 5, et cetera).

poly_roots()

Returns the roots of a polynomial equation, regardless of degree. Unlike the other root-finding functions, it will check for coefficients of zero for the highest power, and 'step down' the degree of the polynomial to the appropriate case. Additionally, it will check for coefficients of zero for the lowest power terms, and add zeros to its root list before calling one of the root-finding functions.

By default, poly_roots() will use the Hessenberg matrix method for solving polynomials. This can be changed by calling "poly_options()".

get_hessenberg() DEPRECATED

Returns 1 or 0 depending upon whether poly_roots() always makes use of the Hessenberg matrix method or not.

NOTE: this function is replaced by the option function poly_option().

set_hessenberg() DEPRECATED

Sets or removes the condition that forces the use of the Hessenberg matrix regardless of the polynomial's degree. A zero argument forces the use of classical methods for polynomials of degree less than five, a non-zero argument forces poly_roots() to always use the matrix method. The default state of the module is to always use the matrix method. This is a complete change from the default behavior in versions less than v2.50.

NOTE: this function is replaced by the option function poly_option().

Classical Functions

These are the functions that solve polynomials via the classical methods. Quartic, cubic, quadratic, and even linear equations may be solved with these functions. They are all exported under the tag "classical".

"poly_roots()" will use these functions if the hessenberg option is set to 0, and if the degree of the polynomial is four or less.

The leading coefficient $a must always be non-zero for the classical functions.

linear_roots()

Here for completeness's sake more than anything else. Returns the solution for

ax + b = 0

by returning -b/a. This may be in either a scalar or an array context.

quadratic_roots()

Gives the roots of the quadratic equation

ax**2 + bx + c = 0

using the well-known quadratic formula. Returns a two-element list.

cubic_roots()

Gives the roots of the cubic equation

ax**3 + bx**2 + cx + d = 0

by the method described by R. W. D. Nickalls (see the "ACKNOWLEDGMENTS" section below). Returns a three-element list. The first element will always be real. The next two values will either be both real or both complex numbers.

quartic_roots()

Gives the roots of the quartic equation

ax**4 + bx**3 + cx**2 + dx + e = 0

using Ferrari's method (see the "ACKNOWLEDGMENTS" section below). Returns a four-element list. The first two elements will be either both real or both complex. The next two elements will also be alike in type.

Sturm Functions

These are the functions that create and make use of the Sturm sequence. They are all exported under the tag "sturm".

poly_real_root_count()

Return the number of unique, real roots of the polynomial.

$unique_roots = poly_real_root_count(@coefficients);

For example, the equation (x + 3)**3 forms the polynomial x**3 + 9x**2 + 27x + 27, but poly_real_root_count(1, 9, 27, 27) will return 1 since all three roots are identical.

Likewise, poly_real_root_count(1, -8, 25) will return 0 because the two roots of x**2 -8x + 25 are both complex.

This function is the all-in-one function to use instead of

my @chain = poly_sturm_chain(@coefficients);

return sturm_sign_count(sturm_sign_minus_inf(\@chain)) -
        sturm_sign_count(sturm_sign_plus_inf(\@chain));

if you don't intend to use the Sturm chain for anything else.

poly_sturm_chain()

Returns the chain of Sturm functions used to evaluate the number of roots of a polynomial in a range of X values.

If you feed in a sequence of X values to the Sturm functions, you can tell where the (real, not complex) roots of the polynomial are by counting the number of times the Y values change sign.

See "poly_real_root_count()" above for an example of its use.

sturm_real_root_range_count()

Return the number of unique, real roots of the polynomial.

my($x0, $x1) = (0, 1000);

my @chain = poly_sturm_chain(@coefficients);
$unique_roots = sturm_real_root_range_count(\@chain, $x0, $x1);

Internally, this is equivalent to:

my @signs = sturm_sign_chain(\@chain, [$x0, $x1]);
return sturm_sign_count(@{$signs[0]}) - sturm_sign_count(@{$signs[1]});

Sturm Sign Sequence Functions

sturm_sign_chain()

sturm_sign_minus_inf()

sturm_sign_plus_inf()

These functions return the array of signs that are used by the functions "poly_real_root_count()" and "sturm_real_root_range_count()" to find the number of real roots in a polynomial.

In normal use you will probably never need to use them, unless you want to examine the internals of the Sturm functions:

#
# Examine the sign changes that occur at each endpoint of
# the x range.
#
my(@coefficients) = (1, 4, 7, 23);
my(@xvals) = (-12, 12);

my @chain = poly_sturm_chain( @coefficients);
my @signs = sturm_sign_chain(\@chain, \@xvals);  # An array of arrays.

print "\nPolynomial: [", join(", ", @coefficients), "]\n";

foreach my $j (0..$#signs)
{
    my @s = @{$signs[$j]};
    print $xval[$j], "\n",
          "\t", join(", ", @s), "], sign count = ",
          sturm_sign_count(@s), "\n\n";
}

Similar examinations can be made at plus and minus infinity:

#
# Examine the sign changes that occur between plus and minus
# infinity.
#
my @coefficients = (1, 4, 7, 23);

my @chain = poly_sturm_chain( @coefficients);
my @smi = sturm_sign_minus_inf(\@chain);
my @spi = sturm_sign_plus_inf(\@chain);

print "\nPolynomial: [", join(", ", @coefficients), "]\n";

print "Minus Inf:\n",
      "\t", join(", ", @smi), "], sign count = ",
      sturm_sign_count(@smi), "\n\n";

print "Plus Inf:\n",
      "\t", join(", ", @spi), "], sign count = ",
      sturm_sign_count(@spi), "\n\n";

sturm_sign_count()

Counts and returns the number of sign changes in a sequence of signs, such as those returned by the "Sturm Sign Sequence Functions"

See "poly_real_root_count()" and "sturm_real_root_range_count()" for examples of its use.

sturm_bisection_roots()

Return the real roots counted by "sturm_real_root_range_count()". Uses the bisection method combined with sturm_real_range_count() to narrow the range to a single root, then uses "laguerre()" to find the value.

my($from, $to) = (-1000, 0);
my @chain = poly_sturm_chain(@coefficients);
my @roots = sturm_bisection_roots(\@chain, $from, $to);

As it is using the Sturm functions, it will find only the real roots.

Utility Functions

These are internal functions used by the other functions listed above, but which may also be useful to the user. They are all exported under the tag "utility".

epsilon()

Returns the machine epsilon value that was calculated when this module was loaded.

The value may be changed, although this in general is not recommended.

my $old_epsilon = epsilon($new_epsilon);

The previous value of epsilon may be saved to be restored later.

The Wikipedia article at http://en.wikipedia.org/wiki/Machine_epsilon/ has more information on the subject.

fltcmp()

Compare two floating point numbers within a degree of accuracy.

Like most functions ending in "cmp", this one returns -1 if the first argument tests as less than the second argument, 1 if the first tests greater than the second, and 0 otherwise. Comparisons are made within a tolerance range that may be set with "poly_tolerance()".

#
# Set a very forgiving comparison tolerance.
#
poly_tolerance(fltcmp => 1e-5);
my @x = poly_roots(@cubic);
my @y = poly_evaluate(\@cubic, \@x);

if (fltcmp($y[0], 0.0) == 0 and
    fltcmp($y[1], 0.0) == 0 and
    fltcmp($y[2], 0.0) == 0)
{
    print "Roots found: (", join(", ", @x), ")\n";
}
else
{
    print "Problem root-finding for [", join(", ", @cubic), "]\n";
}

laguerre()

A numerical method for finding a root of an equation, especially made for polynomials.

@roots = laguerre(\@coefficients, \@xvalues);
push @roots, laguerre(\@coefficients, $another_xvalue);

For each x value the function will attempt to find a root closest to it.

poly_iteration()

Sets the limit to the number of iterations that a solving method may go through before giving up trying to find a root. Each method of root-finding used by "poly_roots()", "sturm_bisection_roots()", and "laguerre()" has its own iteration limit, which may be found, like "poly_option()", by simply looking at the return value of poly_iteration().

#
# Get all of the current iteration limits.
#
my %its_limits = poly_iteration();

#
# Double the limit for the hessenberg method, but set the limit
# for Laguerre's method to 20.
#
poly_iteration(hessenberg => $its_limits{hessenberg} * 2, laguerre => 12);

#
# Reset the limits with the former values, but save the values we had
# for later.
#
my %hl_limits = poly_iteration(%its_limits);

There are iteration limit values for:

hessenberg

The numeric method used by poly_roots(), if the hessenberg option is set. Its default value is 60.

laguerre

The numeric method used by laguerre(). Laguerre's method is used within sturm_bisection_roots() once an individual root has been found within a range, and of course it may be called independently. Its default value is 60.

sturm_bisection

The bisection method used to find roots within a range. Its default value is 100.

poly_tolerance()

Set the degree of accuracy needed for comparisons to be equal or roots to be found. Amongst the root finding functions this currently this only affects laguerre(), as the Hessenberg matrix method determines how close it needs to get using a complicated formula based on epsilon.

#
# Quadruple the tolerance for Laguerre's method.
#
my %tolerances = poly_tolerance();

poly_tolerance(laguerre => 4 * $tolerances{laguerre});

There are iteration limit values for:

laguerre

The numeric method used by laguerre(). Laguerre's method is used within sturm_bisection_roots() once an individual root has been found within a range, and of course it may be called independently. Its default value is 60.

fltcmp

A comparison function that determines if one argument is less than, equal to, or greater than, the other. Comparisons are made within a range determined by the tolerance.

poly_derivative()

@derivative = poly_derivative(@coefficients);

Returns the coefficients of the first derivative of the polynomial. Leading zeros are removed before returning the derivative, so the length of the returned polynomial may be even shorter than expected from the length of the original polynomial. Returns an empty list if the polynomial is a simple constant.

poly_antiderivative()

Returns the coefficients of the antiderivative of the polynomial. The constant term is set to zero; to override this use

@integral = poly_antiderivative(@coefficients);
$integral[$#integral] = $const_term;

simplified_form()

Return the polynomial adjusted by removing any leading zero coefficients and placing it in a monic polynomial form (all coefficients divided by the coefficient of the highest power).

poly_evaluate()

Returns the values of the polynomial given a list of arguments. Unlike most of the above functions, this takes the reference of the coefficient list, which lets the function take a single x-value or multiple x-values passed in as a reference.

The function may return a list...

my @coefficients = (1, -12, 0, 8, 13);
my @xvals = (0, 1, 2, 3, 5, 7);
my @yvals = poly_evaluate(\@coefficients, \@xvals);

print "Polynomial: [", join(", ", @coefficients), "]\n";

for my $j (0..$#yvals) {
    print "Evaluates at ", $xvals[$j], " to ", $yvals[$j], "\n";
}

or return a scalar.

my $x_median = ($xvals[0] + $xvals[$#xvals])/2.0;
my $y_median = poly_evaluate(\@coefficients, $x_median);

poly_nonzero_term_count()

Returns a simple count of the number of coefficients that aren't zero.

poly_constmult()

Simple function to multiply all of the coefficients by a constant. Like poly_evaluate(), uses the reference of the coefficient list.

my @coefficients = (1, 7, 0, 12, 19);
my @coef3 = poly_constmult(\@coefficients, 3);

poly_divide()

Divide one polynomial by another. Like poly_evaluate(), the function takes a reference to the coefficient list. It returns a reference to both a quotient and a remainder.

my @coefficients = (1, -13, 59, -87);
my @polydiv = (3, -26, 59);
my($q, $r) = poly_divide(\@coefficients, \@polydiv);
my @quotient = @$q;
my @remainder = @$r;

EXPORT

There are no default exports. The functions may be individually named in an export list, but there are also four export tags: classical, numeric, sturm, and utility.

ACKNOWLEDGMENTS

The cubic

The cubic is solved by the method described by R. W. D. Nickalls, "A New Approach to solving the cubic: Cardan's solution revealed," The Mathematical Gazette, 77, 354-359, 1993.

Dr. Nickalls was kind enough to send me his article, with notes and revisions, and directed me to a Matlab script that was based on that article, written by Herman Bruyninckx, of the Dept. Mechanical Eng., Div. PMA, Katholieke Universiteit Leuven, Belgium. This function is an almost direct translation of that script, and I owe Herman Bruyninckx for creating it in the first place.

Beginning with version 2.51, Dr. Nikalls's paper is included in the references directory of this package. Dr. Nickalls has also made his paper available at http://www.nickalls.org/dick/papers/maths/cubic1993.pdf.

This article is also available on http://www.2dcurves.com/cubic/cubic.html, and there is a nice discussion of his method at http://www.sosmath.com/algebra/factor/fac111/fac111.html.

Dick Nickalls, dick@nickalls.org

Herman Bruyninckx, Herman.Bruyninckx@mech.kuleuven.ac.be, has web page at http://www.mech.kuleuven.ac.be/~bruyninc. His matlab cubic solver is at http://people.mech.kuleuven.ac.be/~bruyninc/matlab/cubic.m.

Andy Stein has written a version of cubic.m that will work with vectors. It is included with this package in the eg directory.

The quartic

The method for quartic solution is Ferrari's, as described in the web page Karl's Calculus Tutor at http://www.karlscalculus.org/quartic.html. I also made use of some short cuts mentioned in web page Ask Dr. Math FAQ, at http://forum.swarthmore.edu/dr.math/faq/faq.cubic.equations.html.

Quintic and higher.

Back when this module could only solve polynomials of degrees 1 through 4, Matz Kindahl, the author of Math::Polynomial, suggested the poly_roots() function. Later on, Nick Ing-Simmons, who was working on a perl binding to the GNU Scientific Library, sent a perl translation of Hiroshi Murakami's Fortran implementation of the QR Hessenberg algorithm, and it fit very well into the poly_roots() function. Quintics and higher degree polynomials can now be solved, albeit through numeric analysis methods.

Hiroshi Murakami's Fortran routines were at http://netlib.bell-labs.com/netlib/, but do not seem to be available from that source anymore.

He referenced the following articles:

R. S. Martin, G. Peters and J. H. Wilkinson, "The QR Algorithm for Real Hessenberg Matrices", Numer. Math. 14, 219-231(1970).

B. N. Parlett and C. Reinsch, "Balancing a Matrix for Calculation of Eigenvalues and Eigenvectors", Numer. Math. 13, 293-304(1969).

Alan Edelman and H. Murakami, "Polynomial Roots from Companion Matrix Eigenvalues", Math. Comp., v64,#210, pp.763-776(1995).

For starting out, you may want to read

William Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling Numerical Recipes in C. Cambridge University Press, 1988. They have a web site for their book, http://www.nr.com/.

Sturm's Sequence and Laguerre's Method

Dörrie, Heinrich. 100 Great Problems of Elementary Mathematics; Their History and Solution. New York: Dover Publications, 1965. Translated by David Antin.

    Discusses Charles Sturm's 1829 paper with an eye towards mathematical proof rather than an algorithm, but is still very useful.

Glassner, Andrew S. Graphics Gems. Boston: Academic Press, 1990.

    The chapter "Using Sturm Sequences to Bracket Real Roots of Polynomial Equations" (by D. G. Hook and P. R. McAree) has a clearer description of the actual steps needed to implement Sturm's method.

Acton, Forman S. Numerical Methods That Work. New York: Harper & Row, Publishers, 1970.

    Lively, opinionated book on numerical equation solving. I looked it up when it became obvious that everone was quoting Acton when discussing Laguerre's method.

SEE ALSO

Forsythe, George E., Michael A. Malcolm, and Cleve B. Moler Computer Methods for Mathematical Computations. Prentice-Hall, 1977.

AUTHOR

John M. Gamble may be found at jgamble@cpan.org

1 POD Error

The following errors were encountered while parsing the POD:

Around line 2320:

Non-ASCII character seen before =encoding in 'Dörrie,'. Assuming CP1252