Math::Polynomial::Solve - Find the roots of polynomial equations.
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. # set_hessenberg(0); my @x = poly_roots(@coefficients);
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);
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);
use Math::Polynomial::Solve qw(:sturm); # Find the number of unique real roots of the polynomial. my $unique_roots = poly_real_root_count(@coefficients);
This package supplies a set of functions that find the roots of polynomials, along with some utility functions. Polynomials up to the quartic may be solved directly by numerical formulae. Polynomials of fifth and higher powers will be solved by an iterative method, as there are no general solutions for fifth and higher powers.
The linear, quadratic, cubic, and quartic *_roots() functions all expect to have a non-zero value for the $a term.
If the constant term is zero then the first value returned in the list of answers will always be zero, for all functions.
In addition to the root-finding functions, the internal functions have also been exported for your use. See the EXPORT section for a list of functions exported via the :utiltiy tag.
Returns 1 or 0 depending upon whether poly_roots() always makes use of the Hessenberg matrix method or not.
NOTE: this function will probably be deprecated in future releases in favor of a hash-like option function.
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 version 2.50.
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.
If set_substitution() has been called with a non-zero argument, the simplification techniques described in the function description will also be used. With these methods, it is possible to solve a polynomial of degree higher than 4 without using the matrix method, as long as it meets these rather specialized conditions and if there has been a call to set_hessenberg(0).
By default, poly_roots() will use the Hessenberg matrix method for solving polynomials.
If the function set_hessenberg() is called with an argument of 0, poly_roots() attempts to use one of the classical root-finding functions listed below, if the degree of the polynomial is four or less.
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.
-b/a
Gives the roots of the quadratic equation
ax**2 + bx + c = 0
using the well-known quadratic formula. Returns a two-element list.
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.
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.
Return the number of unique, real roots of the polynomial.
$unique_roots = poly_real_root_count(@coefficients);
For example, the polynomial (x + 3)**3 forms the polynomial x**3 + 9x**2 + 27x + 27, but poly_real_root_count() will return 1 for that polynomial since all three roots are the same.
@derivative = poly_derivative(@coefficients);
Returns the coefficients of the first derivative of the polynomial. Because leading zeros are removed before returning the derivative, the resulting polynomial may be reduced by more than just one than the length of the original polynomial. Returns an empty list if the polynomial is a simple constant.
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;
Returns the machine epsilon value that was calculated when this module was loaded. The Wikipedia article at http://en.wikipedia.org/wiki/Machine_epsilon/ has more information on the subject.
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).
Returns the values of the polynomial given a list of arguments. Unlike the above functions, this takes the reference of the coefficient list, in order to feed the multiple x-values (that are also passed in as a reference).
my @polynomial = (1, -12, 0, 8, 13); my @xvals = (0, 1, 2, 3, 5, 7); my @yvals = poly_evaluate(\@polynomial, \@xvals); print "Polynomial: [", join(", ", @polynomial), "]\n"; for my $j (0..$#yvals) { print "Evaluates at ", $xvals[$j], " to ", $yvals[$j], "\n"; }
Simple function to multiply all of the coefficients by a constant. Like poly_evaluate(), uses the reference of the coefficient list.
my @polynomial = (1, 7, 0, 12, 19); my @polynomial3 = poly_constmult(\@polynomial, 3);
Divide one polynomial by another, returning both a quotient and a remainder.
my @polynomial = (1, -13, 59, -87); my @polydiv = (3, -26, 59); my($q, $r) = poly_divide(\@polynomial, \@polydiv); my @quotient = @$q; my @remainder = @$r;
Like poly_evaluate(), the function takes references to the coefficient lists. It returns references to the quotient and the remainder.
There are no default exports. The functions may be named in an export list.
There are also four export tags.
Exports the four root-finding functions for polynomials of degree less than five: linear_roots(), quadratic_roots(), cubic_roots(), quartic_roots().
Exports the root-finding function and the functions that affect its behavior: poly_roots(), get_hessenberg() and set_hessenberg(). set_substitution().
Exports the functions that provide the Sturm functions: poly_sturm_chain() and poly_real_root_count(). Currently only poly_real_root_count() is useful to the typical end user.
Exports the functions that are used internally by other functions, and which might be useful to you as well: epsilon(), poly_evaluate(), poly_derivative(), poly_antiderivatve(), poly_constmult(), poly_divide(), and simplified_form().
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 of this module, Dr. Nikalls's paper is included in the references directory. 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 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.
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
Numerical Recipes in C, by William Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling, Cambridge University Press. They have a web site for their book, http://www.nr.com/.
100 Great Problems of Elementary Mathematics: Their History and Solution, by Heinrich Dorrie, 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.
Graphics Gems, by xxx. 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.
Forsythe, George E., Michael A. Malcolm, and Cleve B. Moler (1977), Computer Methods for Mathematical Computations, Prentice-Hall.
John M. Gamble may be found at jgamble@cpan.org
To install Math::Polynomial::Solve, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Math::Polynomial::Solve
CPAN shell
perl -MCPAN -e shell install Math::Polynomial::Solve
For more information on module installation, please visit the detailed CPAN module installation guide.