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

NAME

Math::Polynom - Operations on polynoms

SYNOPSIS

use Math::Polynom;

To create the polynom '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 polynoms
my $p3 = $p1->multiply(4.5); # multiply a polynom with a constant

my $p3 = $p1->add($p2);      # add 2 polynoms
my $p3 = $p1->add(3.6);      # add a constant to a polynom

my $p3 = $p1->minus($p2);    # substract 2 polynoms
my $p3 = $p1->minus(1.5);    # substract a constant to a polynom

my $p3 = $p1->negate();      # negate a polynom

my $p3 = $p1->divide(3.2);   # divide a polynom by a constant

my $v = $p1->eval(1.35);     # evaluate the polynom on a given value

my $p3 = $p1->derivate();    # return the derivate of a polynom

print $p1->stringify."\n";   # stringify polynom

To try to find a root to a polynom 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 polynoms!! 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 polynom of one variable with up to hundreds of coefficients and non integer powers (ex: powers with decimals). I also needed thorough fault handling. Other CPAN modules providing operations on polynoms did not support those requirements.

If what you need is to manipulate simple polynoms with integer powers, without risks of failures, check out Math::Polynomial since it provides a more complete api than this one.

An instance of Math::Polynom is a representation of a 1-variable polynom. It supports a few basic operations specific to polynoms 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 and Secant 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 polynom.

$p3 = $p1->add($p2)

Return a new polynom that is the sum of the current polynom with the polynom $p2. If $p2 is a scalar, we add it to the current polynom as a numeric constant.

Croaks if provided with weird arguments.

$p3 = $p1->minus($p2)

Return a new polynom that is the current polynom minus the polynom $p2. If $p2 is a scalar, we substract it from the current polynom as a numeric constant.

Croaks if provided with weird arguments.

$p3 = $p1->multiply($p2)

Return a new polynom that is the current polynom multiplied by $p2. If $p2 is a scalar, we multiply all the coefficients in the current polynom with it.

Croaks if provided with weird arguments.

$p3 = $p1->negate()

Return a new polynom in which all coefficients have the negated sign of those in the current polynom.

$p3 = $p1->divide($float)

Return a new polynom in which all coefficients are equal to those of the current polynom divided by the number $float.

Croaks if provided with weird arguments.

$p3 = $p1->derivate()

Return a new polynom that is the derivate of the current polynom.

$v = $p1->eval($float)

Evaluate the current polynom 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 polynom. 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 polynom. Beware that this require your polynom 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 polynom is empty, newton_raphson dies and $p1->error is set to the code Math::Polynom::ERROR_EMPTY_POLYNOM.

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 polynom
            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 polynom. 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 polynom should be continuous. Therefore, the secant method might fail on polynomial 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 polynom is empty, secant dies and $p1->error is set to the code Math::Polynom::ERROR_EMPTY_POLYNOM.

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 polynom. 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 polynom 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 polynome 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 polynom 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.

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 polynom. 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 polynom 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 polynom'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 polynom 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().

ERROR HANDLING

Each method of a polynom 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 polynom'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 polynom 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 polynom (such as newton_raphson)
Math::polynom::ERROR_WRONG_SIGNS means that the polynom 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)

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

SEE ALSO

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

VERSION

$Id: Polynom.pm,v 1.10 2007/01/12 09:23:28 erwan Exp $

THANKS

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

AUTHOR

Erwan Lemonnier <erwan@cpan.org>

COPYRIGHT AND LICENSE

This code is distributed under the same terms as Perl itself.

DISCLAIMER OF WARRANTY

This is free code and comes with no warranty. The author declines any personal responsibility regarding the use of this code or the consequences of its use.