/ 
Math-Polynom-0.10
River stage one • 1 direct dependent • 1 total dependent
/ Math::Polynom

NAME

Math::Polynom - Operations on polynomials

SYNOPSIS

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

DESCRIPTION

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 ughly 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 support 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.

API

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

DEBUG

To display debug information, set in your code:

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

ERROR HANDLING

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

BUGS AND LIMITATIONS

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

REPOSITORY

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 ALSO

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

VERSION

$Id: Polynom.pm,v 1.5 2007/04/17 19:47:27 erwan_lemonnier Exp $

THANKS

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

AUTHORS

Erwan Lemonnier <erwan@cpan.org>, as part of the Pluto developer group at the Swedish Premium Pension Authority.

LICENSE

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.