Erwan Lemonnier >
Math-Polynom >
Math::Polynom

Module Version: 0.13
Math::Polynom - Operations on polynomials

use Math::Polynom;

To create the polynomial 'x^3 + 4*x^2 + 1', write:

my $p1 = Math::Polynom->new(3 => 1, 2 => 4, 0 => 1);

To create '3.5*x^4.2 + 1.78*x^0.9':

my $p2 = Math::Polynom->new(4.2 => 3.5, 0.9 => 1.78);

Common operations:

my $p3 = $p1->multiply($p2); # multiply 2 polynomials my $p3 = $p1->multiply(4.5); # multiply a polynomial with a constant my $p3 = $p1->add($p2); # add 2 polynomials my $p3 = $p1->add(3.6); # add a constant to a polynomial my $p3 = $p1->minus($p2); # substract 2 polynomials my $p3 = $p1->minus(1.5); # substract a constant to a polynomial my $p3 = $p1->negate(); # negate a polynomial my $p3 = $p1->divide(3.2); # divide a polynomial by a constant my $v = $p1->eval(1.35); # evaluate the polynomial on a given value my $p3 = $p1->derivate(); # return the derivate of a polynomial print $p1->stringify."\n"; # stringify polynomial

To try to find a root to a polynomial using the Newton Raphson method:

my $r; eval { $r = $p1->newton_raphson(guess => 2, precision => 0.001); }; if ($@) { if ($p1->error) { # that's an internal error if ($p1->error == Math::Polynom::ERROR_NAN) { # bumped on a complex number } } else { # either invalid arguments (or a bug in solve()) } }

Same with the secant method:

eval { $r = $p1->secant(p0 => 0, p2 => 2, precision => 0.001); };

What! Yet another module to manipulate polynomials!! No, don't worry, there is a good reason for this one ;)

I needed (for my work at a large financial institution) a robust way to compute the internal rate of return (IRR) of various cashflows. An IRR is typically obtained by solving a usually ugly looking polynomial of one variable with up to hundreds of coefficients and non integer powers (ex: powers with decimals). I also needed thorough exception handling. Other CPAN modules providing operations on polynomials did not fill those requirements.

If what you need is to manipulate simple polynomials with integer powers, without concern for failures, check out Math::Polynomial since it provides a more complete api than Math::Polynom.

An instance of Math::Polynom is a representation of a 1-variable polynomial. It supports a few basic operations specific to polynomials such as addition, substraction and multiplication.

Math::Polynom also implements various root finding algorithms (which is kind of the main purpose of this module) such as the Newton Raphson, Secant and Brent methods.

- $p1 =
**new(%power_coef)** -
Create a new Math::Polynom. Each key in the hash

`%power_coef`

is a power and each value the corresponding coefficient. - $p3 = $p1->
**clone()** -
Return a clone of the current polynomial.

- $p3 = $p1->
**add($p2)** -
Return a new polynomial that is the sum of the current polynomial with the polynomial

`$p2`

. If`$p2`

is a scalar, we add it to the current polynomial as a numeric constant.Croaks if provided with weird arguments.

- $p3 = $p1->
**minus($p2)** -
Return a new polynomial that is the current polynomial minus the polynomial

`$p2`

. If`$p2`

is a scalar, we substract it from the current polynomial as a numeric constant.Croaks if provided with weird arguments.

- $p3 = $p1->
**multiply($p2)** -
Return a new polynomial that is the current polynomial multiplied by

`$p2`

. If`$p2`

is a scalar, we multiply all the coefficients in the current polynomial with it.Croaks if provided with weird arguments.

- $p3 = $p1->
**negate()** -
Return a new polynomial in which all coefficients have the negated sign of those in the current polynomial.

- $p3 = $p1->
**divide($float)** -
Return a new polynomial in which all coefficients are equal to those of the current polynomial divided by the number

`$float`

.Croaks if provided with weird arguments.

- $p3 = $p1->
**derivate()** -
Return a new polynomial that is the derivate of the current polynomial.

- $v = $p1->
**eval($float)** -
Evaluate the current polynomial on the value

`$float`

.If you call

`eval`

with a negative value that would yield a complex (non real) result,`eval`

will no complain but return the string 'nan'.Croaks if provided with weird arguments.

- $s = $p1->
**stringify()** -
Return a basic string representation of the current polynomial. For exemple '3*x^5 + 2*x^2 + 1*x^0'.

- $r = $p1->
**newton_raphson(guess => $float1, precision => $float2, max_depth => $integer)** -
Uses the Newton Raphson algorithm to approximate a root for this polynomial. Beware that this require your polynomial AND its derivate to be continuous. Starts the search with

`guess`

and returns the root when the difference between two consecutive estimations of the root is smaller than`precision`

. Make at most`max_depth`

iterations.If

`guess`

is omitted, 1 is used as default. If`precision`

is omitted, 0.1 is used as default. If`max_depth`

is omitted, 100 is used as default.`newton_raphson`

will fail (croak) in a few cases: If the successive approximations of the root still differ with more than`precision`

after`max_depth`

iterations,`newton_raphson`

dies, and`$p1->error`

is set to the code Math::Polynom::ERROR_MAX_DEPTH. If an approximation is not a real number,`newton_raphson`

dies and`$p1->error`

is set to the code Math::Polynom::ERROR_NAN. If the polynomial is empty,`newton_raphson`

dies and`$p1->error`

is set to the code Math::Polynom::ERROR_EMPTY_POLYNOM. If`newton_raphson`

converges toward a number but this number is not a root (ie the polynomial evaluates to a large number on it),`newton_raphson`

dies and the error attribute is set to Math::Polynom::ERROR_NOT_A_ROOT.`newton_raphson`

will also croak if provided with weird arguments.Exemple:

eval { $p->newton_raphson(guess => 1, precision => 0.0000001, max_depth => 50); }; if ($@) { if ($p->error) { if ($p->error == Math::Polynom::ERROR_MAX_DEPTH) { # do something wise } elsif ($p->error == Math::Polynom::ERROR_MAX_DEPTH) { # do something else } else { # empty polynomial die "BUG!"; } } else { die "newton_raphson died for unknown reason"; } }

- $r = $p1->
**secant(p0 => $float1, p1 => $float2, precision => $float3, max_depth => $integer)** -
Use the secant method to approximate a root for this polynomial.

`p0`

and`p1`

are the two start values to initiate the search,`precision`

and`max_depth`

have the same meaning as for`newton_raphson`

.The polynomial should be continuous. Therefore, the secant method might fail on polynomialial having monoms with degrees lesser than 1.

If

`precision`

is omitted, 0.1 is used as default. If`max_depth`

is omitted, 100 is used as default.`secant`

will fail (croak) in a few cases: If the successive approximations of the root still differ with more than`precision`

after`max_depth`

iterations,`secant`

dies, and`$p1->error`

is set to the code Math::Polynom::ERROR_MAX_DEPTH. If an approximation is not a real number,`secant`

dies and`$p1->error`

is set to the code Math::Polynom::ERROR_NAN. If the polynomial is empty,`secant`

dies and`$p1->error`

is set to the code Math::Polynom::ERROR_EMPTY_POLYNOM. If`secant`

converges toward a number but this number is not a root (ie the polynomial evaluates to a large number on it),`secant`

dies and the error attribute is set to Math::Polynom::ERROR_NOT_A_ROOT.`secant`

will also croak if provided with weird arguments. - $r = $p1->
**brent(a => $float1, b => $float2, precision => $float3, max_depth => $integer)** -
Use Brent's method to approximate a root for this polynomial.

`a`

and`b`

are two floats such that`p1->eval(a)`

and`p1->eval(b)`

have opposite signs.`precision`

and`max_depth`

have the same meaning as for`newton_raphson`

.The polynomial should be continuous on the interval [a,b].

Brent's method is considered to be one of the most robust root finding methods. It alternatively uses the secant, inverse quadratic interpolation and bisection to find the next root candidate at each iteration, making it a robust but quite fast converging method.

The difficulty with Brent's method consists in finding the start values a and b for which the polynomial evaluates to opposite signs. This is somewhat simplified in Math::Polynom by the fact that

`eval()`

automatically sets`xpos()`

and`xneg()`

when possible.If

`precision`

is omitted, 0.1 is used as default. If`max_depth`

is omitted, 100 is used as default.`brent`

will fail (croak) in a few cases: If the successive approximations of the root still differ with more than`precision`

after`max_depth`

iterations,`brent`

dies, and`$p1->error`

is set to the code Math::Polynom::ERROR_MAX_DEPTH. If an approximation is not a real number,`brent`

dies and`$p1->error`

is set to the code Math::Polynom::ERROR_NAN. If the polynomial is empty,`brent`

dies and`$p1->error`

is set to the code Math::Polynom::ERROR_EMPTY_POLYNOM. If provided with a and b that does not lead to values having opposite signs,`brent`

dies and`$p1->error`

is set to the code Math::Polynom::ERROR_WRONG_SIGNS. If`brent`

converges toward a number but this number is not a root (ie the polynomial evaluates to a large number on it),`brent`

dies and the error attribute is set to Math::Polynom::ERROR_NOT_A_ROOT.`brent`

will also croak if provided with weird arguments. - $p1->
**error**, $p1->**error_message** -
Respectively the error code and error message set by the last method that failed to run on this polynomial. For exemple, if

`newton_raphson`

died, you would access the code of the error with`error()`

and a message describing the context of the error in details with`error_message`

.If the polynomial has no error,

`error`

returns Math::polynom::NO_ERROR and`error_message`

returns undef. - $p1->
**iterations** -
Return the number of iterations it took to find the polynomial's root. Must be called after calling one of the root finding methods.

- $p1->
**xpos**, $p1->**xneg** -
Each time

`eval`

is called, it checks whether we know a value xpos for which the polynomial evaluates to a positive value. If not and if the value provided to`eval`

lead to a positive result, this value is stored in`xpos`

. Same thing with`xneg`

and negative results.This comes in handy when you wish to try the Brent method after failing with the secant or Newton methods. If you are lucky, those failed attempts will have identified both a xpos and xneg that you can directly use as a and b in

`brent()`

.

To display debug information, set in your code:

local $Math::Polynom::DEBUG = 1;

Each method of a polynomial may croak if provided with wrong arguments. Methods that take arguments do thorough controls on whether the arguments are of the proper type and in the right quantity. If the error is internal, the method will croak after setting the polynomial's error and error_message to specific values.

Math::Polynom defines a few error codes, returned by the method `error`

:

**Math::polynom::NO_ERROR**is the default return value of method`error`

, and is always set to 0.**Math::polynom::ERROR_NAN**means the function jammed on a complex number. Most likely because your polynomial is not continuous on the search interval.**Math::polynom::ERROR_DIVIDE_BY_ZERO**means what it says.**Math::polynom::ERROR_MAX_DEPTH**means the root finding algorithm failed to find a good enough root after the specified maximum number of iterations.**Math::polynom::ERROR_EMPTY_POLYNOM**means you tried to perform an operation on an empty polynomial (such as`newton_raphson)`

**Math::polynom::ERROR_WRONG_SIGNS**means that the polynomial evaluates to values having the same signs instead of opposite signs on the boundaries of the interval you provided to start the search of the root (ex: Brent's method)**Math::polynom::ERROR_NOT_A_ROOT**means the root finding method converged toward one value but this value appears not to be a root. A value is accepted as a root if the polynomial evaluates on it to a number between -1 and 1 (ie close enough to 0).

This module is built for robustness in order to run in requiring production environments. Yet it has one limitation: due to Perl's inability at handling large floats, root finding algorithms will get lost if starting on a guess value that is too far from the root. Example:

my $p = Math::Polynom->new(2 => 1, 1 => -2, 0 => 1); # x^2 -2*x +1 $p->newton_raphson(guess => 100000000000000000); # returns 1e17 as the root

The source of Math::Polynom is hosted at sourceforge as part of the xirr4perl project. You can access it at https://sourceforge.net/projects/xirr4perl/.

See Math::Calculus::NewtonRaphson, Math::Polynomial, Math::Function::Roots.

$Id: Polynom.pm,v 1.10 2007/07/11 13:01:48 erwan_lemonnier Exp $

Thanks to Spencer Ogden who wrote the implementation of the Secant algorithm in his module Math::Function::Roots.

Erwan Lemonnier `<erwan@cpan.org>`

, as part of the Pluto developer group at the Swedish Premium Pension Authority.

This code was developed at the Swedish Premium Pension Authority as part of the Authority's software development activities. This code is distributed under the same terms as Perl itself. We encourage you to help us improving this code by sending feedback and bug reports to the author(s).

This code comes with no warranty. The Swedish Premium Pension Authority and the author(s) decline any responsibility regarding the possible use of this code or any consequence of its use.

syntax highlighting: