Math-Polynom

 view release on metacpan or  search on metacpan

examples/newton_raphson.pl  view on Meta::CPAN

while (abs($new_guess - $old_guess) > $precision) {
    $old_guess = $new_guess;

    my $dividend = $derivate->eval($old_guess);

    die "division by zero: polynomial's derivate is 0 at $old_guess"
	if ($dividend == 0);

    $new_guess = $old_guess - $p->eval($old_guess)/$dividend;

    $p->iterations($p->iterations + 1);

    die "reached maximum number of iterations [$max_depth] without getting close enough to the root."
	if ($p->iterations > $max_depth);
}

print "the root of:\n".$p->stringify."\nis: ".$new_guess."\n";

lib/Math/Polynom.pm  view on Meta::CPAN

package Math::Polynom;

use 5.006;
use strict;
use warnings;
use Carp qw(confess croak);
use Data::Dumper;
use Data::Float qw(float_is_nan);
use Scalar::Util qw(looks_like_number);

use accessors qw(error error_message iterations xpos xneg);

use constant NO_ERROR             => 0;
use constant ERROR_NAN            => 1;
use constant ERROR_MAX_DEPTH      => 2;
use constant ERROR_EMPTY_POLYNOM  => 3;
use constant ERROR_DIVIDE_BY_ZERO => 4;
use constant ERROR_WRONG_SIGNS    => 5;
use constant ERROR_NOT_A_ROOT     => 6;

our $VERSION = '0.13';

lib/Math/Polynom.pm  view on Meta::CPAN

    croak __PACKAGE__." ERROR: $msg\n";
}

sub _exception {
    my ($self,$code,$msg,$args) = @_;

    $msg = "ERROR: $msg\nwith polynom:\n".$self->stringify."\n";
    if (defined $args) {
	$msg .= "with arguments:\n".Dumper($args);
    }
    $msg .= "at iteration ".$self->iterations."\n";

    $self->error_message($msg);
    $self->error($code);

    croak $self->error_message;
}

#################################################################
#
#

lib/Math/Polynom.pm  view on Meta::CPAN


    _error("new() got odd number of arguments. can not be a hash") if (scalar(@args) % 2);

    my %hash = @args;
    foreach my $n (@args) {
	_error("at least one argument of new() is not numeric:\n".Dumper(\%hash)) if (!looks_like_number($n));
    }

    my $self = bless({polynom => \%hash},$pkg)->_clean;
    $self->error(NO_ERROR);
    $self->iterations(0);

    return $self;
}

#----------------------------------------------------------------
#
#   clone - return a clone of self
#

sub clone {

lib/Math/Polynom.pm  view on Meta::CPAN

#
#   newton_raphson - attempt to find a polynomial's root with Newton Raphson
#

sub newton_raphson {
    my($self,%hash) = @_;
    my $new_guess = 1;
    my $precision = 0.1;
    my $max_depth = 100;

    $self->iterations(0);
    $self->error(NO_ERROR);

    $new_guess = $hash{guess}     if (exists $hash{guess});
    $precision = $hash{precision} if (exists $hash{precision});
    $max_depth = $hash{max_depth} if (exists $hash{max_depth});

    _error("newton_raphson() got undefined guess")       if (!defined $new_guess);
    _error("newton_raphson() got undefined precision")   if (!defined $precision);
    _error("newton_raphson() got undefined max_depth")   if (!defined $max_depth);
    _error("newton_raphson() got non numeric guess")     if (!looks_like_number($new_guess));

lib/Math/Polynom.pm  view on Meta::CPAN


    while (abs($new_guess - $old_guess) > $precision) {
	$old_guess = $new_guess;

	my $dividend = $derivate->eval($old_guess);
	$self->_exception(ERROR_DIVIDE_BY_ZERO,"division by zero: polynomial's derivate is 0 at $old_guess",\%hash)
	    if ($dividend == 0);

	$new_guess = $old_guess - $self->eval($old_guess)/$dividend;

	$self->iterations($self->iterations + 1);
	$self->_exception(ERROR_MAX_DEPTH,"reached maximum number of iterations [$max_depth] without getting close enough to the root.",\%hash)
	    if ($self->iterations > $max_depth);

	$self->_exception(ERROR_NAN,"new guess is not a real number in newton_raphson().",\%hash)
	    if (float_is_nan($new_guess));
    }

    if (!$self->_is_root($new_guess)) {
	$self->_exception(ERROR_NOT_A_ROOT,"newton_raphson() converges toward $new_guess but that doesn't appear to be a root.",\%hash);
    }

    return $new_guess;

lib/Math/Polynom.pm  view on Meta::CPAN

#
#   secant - implement the Secant algorithm to approximate the root of this polynomial
#

sub secant {
    my ($self,%hash) = @_;
    my $precision = 0.1;
    my $max_depth = 100;
    my ($p0,$p1);

    $self->iterations(0);
    $self->error(NO_ERROR);

    $precision = $hash{precision} if (exists $hash{precision});
    $max_depth = $hash{max_depth} if (exists $hash{max_depth});
    $p0        = $hash{p0}        if (exists $hash{p0});
    $p1        = $hash{p1}        if (exists $hash{p1});

    _error("secant() got undefined precision")      if (!defined $precision);
    _error("secant() got undefined max_depth")      if (!defined $max_depth);
    _error("secant() got non numeric precision")    if (!looks_like_number($precision));

lib/Math/Polynom.pm  view on Meta::CPAN

    my $p;

    $self->_exception(ERROR_NAN,"q0 or q1 are not a real number in first eval in secant()",\%hash)
	if (float_is_nan($q0) || float_is_nan($q1));

    return $p0 if ($q0 == 0);
    return $p1 if ($q1 == 0);

    for (my $depth = 1; $depth <= $max_depth; $depth++) {

	$self->iterations($depth);

	$self->_exception(ERROR_DIVIDE_BY_ZERO,"division by zero with p0=$p0, p1=$p1, q1=q0=$q1 in secant()",\%hash)
	    if (($q1 - $q0) == 0);

	$p = ($q1 * $p0 - $p1 * $q0) / ($q1 - $q0);

	$self->_exception(ERROR_NAN,"p is not a real number in secant()",\%hash)
	    if (float_is_nan($p));

	my $debug = "secant at depth ".$self->iterations.", p0=$p0, p1=$p1, p=$p";

	$p0 = $p1;
	$q0 = $q1;
	$q1 = $self->eval($p);

	$self->_exception(ERROR_NAN,"q1 is not a real number in secant()",\%hash)
	    if (float_is_nan($q1));

	_debug($debug.", poly(p)=$q1");

	if ($q1 == 0 || abs($p - $p1) <= $precision) {
	    if (!$self->_is_root($p)) {
		$self->_exception(ERROR_NOT_A_ROOT,"secant() converges toward $p but that doesn't appear to be a root.",\%hash);
	    }
	    return $p;
	}

	$p1 = $p;
    }

    $self->_exception(ERROR_MAX_DEPTH,"reached maximum number of iterations [$max_depth] without getting close enough to the root in secant()",\%hash);
}

#----------------------------------------------------------------
#
#   brent - implement Brent's method to approximate the root of this polynomial
#

sub brent {
    my($self,%hash) = @_;
    my $precision = 0.1;
    my $max_depth = 100;
    my $mflag;
    my($a,$b,$c,$s,$d);
    my($f_a,$f_b,$f_c,$f_s);

    $self->iterations(0);
    $self->error(NO_ERROR);

    $precision = $hash{precision} if (exists $hash{precision});
    $max_depth = $hash{max_depth} if (exists $hash{max_depth});
    $a         = $hash{a}         if (exists $hash{a});
    $b         = $hash{b}         if (exists $hash{b});

    _error("brent() got undefined precision")      if (!defined $precision);
    _error("brent() got undefined max_depth")      if (!defined $max_depth);
    _error("brent() got non numeric precision")    if (!looks_like_number($precision));

lib/Math/Polynom.pm  view on Meta::CPAN

    }

    $c = $a;
    $f_c = $f_a;

    $mflag = 1;

    # repeat while we haven't found the root nor are close enough to it
    while ($f_b != 0 && abs($b - $a) > $precision) {

	# did we reach the maximum number of iterations?
	$self->_exception(ERROR_MAX_DEPTH,"reached maximum number of iterations [$max_depth] without getting close enough to the root in brent()",\%hash)
	    if ($self->iterations > $max_depth);

	# evaluate f(a), f(b) and f(c) if necessary
	if ($self->iterations != 0) {
	    $f_a = $self->eval($a);
	    $f_b = $self->eval($b);
	    $f_c = $self->eval($c);

	    $self->_exception(ERROR_NAN,"polynomial leads to an imaginary number on a=$a in brent()",\%hash) if (float_is_nan($f_a));
	    $self->_exception(ERROR_NAN,"polynomial leads to an imaginary number on b=$b in brent()",\%hash) if (float_is_nan($f_b));
	    $self->_exception(ERROR_NAN,"polynomial leads to an imaginary number on c=$c in brent()",\%hash) if (float_is_nan($f_c));
	}

	my $debug = "brent at depth ".$self->iterations.", a=$a, b=$b";

	# calculate the next root candidate
	if ($f_a == $f_b) {
	    # we should not be able to get $f_b == $f_a since it's a prerequisite of the method. that would be a bug
	    _error("BUG: got same values for polynomial at a=$a and b=$b:\n".$self->stringify);

	} elsif ( ($f_a != $f_c) && ($f_b != $f_c) ) {
	    # use quadratic interpolation
	    $s = ($a*$f_b*$f_c)/(($f_a - $f_b)*($f_a - $f_c)) +
		($b*$f_a*$f_c)/(($f_b - $f_a)*($f_b - $f_c)) +

lib/Math/Polynom.pm  view on Meta::CPAN

	    $f_a = $f_s;
	}

	# eventually swap $a and $b
	if (abs($f_a) < abs($f_b)) {
	    # in the special case when
	    ($a,$b) = ($b,$a);
	    ($f_a,$f_b) = ($f_b,$f_a);
	}

	$self->iterations($self->iterations + 1);
    }

    if (!$self->_is_root($b)) {
	$self->_exception(ERROR_NOT_A_ROOT,"brent() converges toward $b but that doesn't appear to be a root.",\%hash);
    }

    return $b;
}

1;

lib/Math/Polynom.pm  view on Meta::CPAN

=item $s = $p1->B<stringify()>

Return a basic string representation of the current polynomial. For exemple '3*x^5 + 2*x^2 + 1*x^0'.

=item $r = $p1->B<< 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 C<guess> and returns the root when the difference between two
consecutive estimations of the root is smaller than C<precision>. Make at most C<max_depth>
iterations.

If C<guess> is omitted, 1 is used as default.
If C<precision> is omitted, 0.1 is used as default.
If C<max_depth> is omitted, 100 is used as default.

C<newton_raphson> will fail (croak) in a few cases: If the successive approximations of the root
still differ with more than C<precision> after C<max_depth> iterations, C<newton_raphson> dies,
and C<< $p1->error >> is set to the code Math::Polynom::ERROR_MAX_DEPTH. If an approximation
is not a real number, C<newton_raphson> dies and C<< $p1->error >> is set to the code Math::Polynom::ERROR_NAN.
If the polynomial is empty, C<newton_raphson> dies and C<< $p1->error >> is set to the code
Math::Polynom::ERROR_EMPTY_POLYNOM. If C<newton_raphson> converges toward a number but this number is
not a root (ie the polynomial evaluates to a large number on it), C<newton_raphson> dies and the
error attribute is set to Math::Polynom::ERROR_NOT_A_ROOT.

C<newton_raphson> will also croak if provided with weird arguments.

Exemple:

lib/Math/Polynom.pm  view on Meta::CPAN

Use the secant method to approximate a root for this polynomial. C<p0> and C<p1> are the two start values
to initiate the search, C<precision> and C<max_depth> have the same meaning as for C<newton_raphson>.

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

If C<precision> is omitted, 0.1 is used as default.
If C<max_depth> is omitted, 100 is used as default.

C<secant> will fail (croak) in a few cases: If the successive approximations of the root
still differ with more than C<precision> after C<max_depth> iterations, C<secant> dies,
and C<< $p1->error >> is set to the code Math::Polynom::ERROR_MAX_DEPTH. If an approximation
is not a real number, C<secant> dies and C<< $p1->error >> is set to the code Math::Polynom::ERROR_NAN.
If the polynomial is empty, C<secant> dies and C<< $p1->error >> is set to the code
Math::Polynom::ERROR_EMPTY_POLYNOM. If C<secant> converges toward a number but this number is
not a root (ie the polynomial evaluates to a large number on it), C<secant> dies and the
error attribute is set to Math::Polynom::ERROR_NOT_A_ROOT.

C<secant> will also croak if provided with weird arguments.


lib/Math/Polynom.pm  view on Meta::CPAN

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 C<eval()> automatically sets C<xpos()> and C<xneg()> when possible.

If C<precision> is omitted, 0.1 is used as default.
If C<max_depth> is omitted, 100 is used as default.

C<brent> will fail (croak) in a few cases: If the successive approximations of the root
still differ with more than C<precision> after C<max_depth> iterations, C<brent> dies,
and C<< $p1->error >> is set to the code Math::Polynom::ERROR_MAX_DEPTH. If an approximation
is not a real number, C<brent> dies and C<< $p1->error >> is set to the code Math::Polynom::ERROR_NAN.
If the polynomial is empty, C<brent> dies and C<< $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, C<brent> dies and C<< $p1->error >> is set to the code Math::Polynom::ERROR_WRONG_SIGNS.
If C<brent> converges toward a number but this number is not a root (ie the polynomial evaluates
to a large number on it), C<brent> dies and the
error attribute is set to Math::Polynom::ERROR_NOT_A_ROOT.

C<brent> will also croak if provided with weird arguments.

lib/Math/Polynom.pm  view on Meta::CPAN


Respectively the error code and error message set by the last method that failed to run
on this polynomial. For exemple, if C<newton_raphson> died, you would access the code of the error
with C<error()> and a message describing the context of the error in details with
C<error_message>.

If the polynomial has no error, C<error> returns Math::polynom::NO_ERROR and
C<error_message> returns undef.


=item $p1->B<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.


=item $p1->B<xpos>, $p1->B<xneg>

Each time C<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 C<eval> lead to a positive
result, this value is stored in C<xpos>. Same thing with C<xneg> and negative results.

This comes in handy when you wish to try the Brent method after failing with the secant

lib/Math/Polynom.pm  view on Meta::CPAN

Math::Polynom defines a few error codes, returned by the method C<error>:

=over 4

=item B<Math::polynom::NO_ERROR> is the default return value of method C<error>, and is always set to 0.

=item B<Math::polynom::ERROR_NAN> means the function jammed on a complex number. Most likely because your polynomial is not continuous on the search interval.

=item B<Math::polynom::ERROR_DIVIDE_BY_ZERO> means what it says.

=item B<Math::polynom::ERROR_MAX_DEPTH> means the root finding algorithm failed to find a good enough root after the specified maximum number of iterations.

=item B<Math::polynom::ERROR_EMPTY_POLYNOM> means you tried to perform an operation on an empty polynomial (such as C<newton_raphson)>

=item B<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)

=item B<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 t...

=back

=head1 BUGS AND LIMITATIONS

t/02_test_new.t  view on Meta::CPAN


my $p1;

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

is(ref $p1, "Math::Polynom", "check type");
is_deeply($p1,
	  {
	      polynom => { 1 => 2, 3 => 4 },
	      -error  => 0,
	      -iterations => 0,
	  },
	  "check content"
	);

$p1 = new Math::Polynom(1.8 => 5.2, 3.1 => 0);

is(ref $p1, "Math::Polynom", "check type");
is_deeply($p1,
	  {
	      polynom => { 1.8 => 5.2 },
	      -error  => 0,
	      -iterations => 0,
	  },
	  "check content"
	);

$p1 = new Math::Polynom();

is(ref $p1, "Math::Polynom", "check type (empty polynom)");
is_deeply($p1,{ polynom => {}, -error => 0, -iterations => 0 },"check content (empty polynom)");
ok(!$p1->error,"polynom contains no error");

# test _is_number, indirectly
eval { 
    Math::Polynom->new(1 => 2.5,
		       -1.234728347568237456 => +1.345e23,
		       2 => -1.345e-23,
		       3 => -1.345E-23,
		       4 => -1.345E+23,
		       );

t/12_test_newton_raphson.t  view on Meta::CPAN


sub test_newton_raphson {
    my($p,$args,$want)=@_;
    my $precision = $args->{precision} || 0.1;
    my $v = $p->newton_raphson(%$args);
    ok(alike($v,$want,$precision), $p->stringify." ->newton_raphson(".( (exists $args->{guess}) ? "guess => ".$args->{guess}.", ":"" )."precision => $precision) = $want (got $v)");
}

my $p1 = Math::Polynom->new(2 => 1, 0 => -4);
test_newton_raphson($p1, {}, 2);
is( $p1->iterations, 3, "p1->iterations=3 after search" );
test_newton_raphson($p1, {guess => 1}, 2);
test_newton_raphson($p1, {guess => -1}, -2);
test_newton_raphson($p1, {guess => 100}, 2);
test_newton_raphson($p1, {guess => 100, precision => 0.0001}, 2);
test_newton_raphson($p1, {guess => 100, precision => 0.0000001}, 2);
test_newton_raphson($p1, {guess => -10000, precision => 0.0000001}, -2);

my $p2 = Math::Polynom->new(5 => 5, 3.2 => 4, 0.9 => -2);  # 5*x^5 + 4*x^3.2 - 2*x^0.9
test_newton_raphson($p2, {precision => 0.000000000000001}, 0.6161718040343);

t/12_test_newton_raphson.t  view on Meta::CPAN

test_newton_raphson($p2, {guess => 50, precision => 0.01}, 0.6161718040343);

$v = undef;
eval { $v = $p2->newton_raphson(guess => 1, precision => 0.000000000000001, max_depth => 2) };

if (defined $v) {
    ok(looks_like_number($v) && abs($v - 0.6161718040343) < 0.000000000001, "\$p2 with limited max_depth was solved on this platform [\$v = $v]");
    ok(1,"nop");
    ok(1,"nop");
} else {
    ok( defined $@ && $@ =~ /reached maximum number of iterations/, "newton_raphson() fails on \$p2 with limited max_depth" );
    ok( defined $p2->error_message && $p2->error_message =~ /reached maximum number of iterations/, "\$p2->error_message looks good" );
    is( $p2->error, Math::Polynom::ERROR_MAX_DEPTH, "\$p2->error looks good" );
}

my $p3 = Math::Polynom->new(2 => 1, 1 => -2, 0 => 1); # x^2 -2*x +1
test_newton_raphson($p3, {guess => -10}, 1);
test_newton_raphson($p3, {guess => 10}, 1);
test_newton_raphson($p3, {guess => 1000000}, 1);

# TODO: handle calculation overflows...
#my $v = $p3->newton_raphson(guess => 100000000000000000);

t/14_test_secant.t  view on Meta::CPAN

}

sub test_secant {
    my ($p,$args,$want) = @_;
    my $precision = $args->{precision} || 0.1;
    my $v = $p->secant(%$args);
    ok(alike($v,$want,$precision), $p->stringify." ->secant(p0 => ".$args->{p0}.", p1 => ".$args->{p1}.", precision => $precision) = $want (got $v)");
}

my $p1 = Math::Polynom->new(2 => 1, 0 => -4);
is( $p1->iterations, 0, "p1->iterations=0 before search" );
test_secant($p1, {p0 => 0.5, p1 => 1}, 2);
is( $p1->iterations, 4, "p1->iterations=3 after search" );
test_secant($p1, {p0 => 0.5, p1 => 1}, 2);
test_secant($p1, {p0 => -5, p1 => -1}, -2);
test_secant($p1, {p0 => 0.5, p1 => 1}, 2);
test_secant($p1, {p0 => 0.5, p1 => 1, precision => 0.0001}, 2);
test_secant($p1, {p0 => 0.5, p1 => 1, precision => 0.0000001}, 2);
test_secant($p1, {p0 => 0, p1 => -100, precision => 0.0000001}, -2);
is( $p1->iterations, 15, "p1->iterations=14 after search" );
test_secant($p1, {p0 => 0.5, p1 => 1}, 2);
is( $p1->iterations, 4, "p1->iterations=3 after simpler search" );


my $p2 = Math::Polynom->new(5 => 5, 3.2 => 4, 0.9 => -2);  # 5*x^5 + 4*x^3.2 - 2*x^0.9
test_secant($p2, {p0 => 0.5, p1 => 1, precision => 0.000000000000001}, 0.6161718040343);

eval { test_secant($p2, {p0 => 0.5, p1 => 1, precision => 0.000000000000001}, 0.6161718040343); };
ok( !defined $@ || $@ eq '', "secant() does not fails on polynom 2 with negative guess (newton_raphson would)" );

my $p3 = Math::Polynom->new(2 => 1, 1 => -2, 0 => 1); # x^2 -2*x +1
test_secant($p3, {p0 => 0.5, p1 => 1}, 1);
test_secant($p3, {p0 => 500, p1 => -500}, 1);
test_secant($p3, {p0 => 100000, p1 => 99999}, 1);

# TODO: handle calculation overflows...
my $v;
eval { $v = $p3->secant(p0 => 100000000000000000, p1 => 999999999999999999999, max_depth => 1); };
ok( defined $@ && $@ =~ /reached maximum number of iterations/, "secant() fails when max_depth reached" .((defined $v)?" (v=$v)":"") );
ok( defined $p3->error_message && $p3->error_message =~ /reached maximum number of iterations/, "\$p3->error_message looks good" );
is( $p3->error, Math::Polynom::ERROR_MAX_DEPTH, "\$p3->error looks good" );

# empty polynom error
my $p4 = Math::Polynom->new();
$v = undef;
eval { $v = $p4->secant(p0 => 0, p1 => 1); };
ok( defined $@ && $@ =~ /empty polynom/, "secant() fails on empty polynom".((defined $v)?" (v=$v)":"") );
ok( defined $p4->error_message && $p4->error_message =~ /empty polynom/, "\$p4->error_message looks good" );
is( $p4->error, Math::Polynom::ERROR_EMPTY_POLYNOM, "\$p4->error looks good" );

t/15_test_brent.t  view on Meta::CPAN

sub test_brent {
    my($p,$args,$want)=@_;
    my $precision = $args->{precision} || 0.1;
    my $sym = $args->{sym} || 0;
    my $v = $p->brent(%$args);
    ok(alike($v,$want,$precision,$sym), $p->stringify." ->brent(a => ".$args->{a}.", b => ".$args->{b}.", precision => $precision) = $want (got $v)");
}

# the exemple on wikipedia:
my $p0 = Math::Polynom->new(3 => 1, 2 => 1, 1 => -5, 0 => 3);
is($p0->iterations,0,"p0->iterations=0");
test_brent($p0, {a => -4, b => 4/3}, -3);
is($p0->iterations,8,"p0->iterations=8 after search");

# p1 is 'x^2 - 4'
my $p1 = Math::Polynom->new(2 => 1, 0 => -4);
test_brent($p1, {a => 1, b => 3}, 2);
test_brent($p1, {a => 3, b => 1}, 2);
test_brent($p1, {a => -5, b => -1}, -2);
test_brent($p1, {a => -1, b => -5}, -2);
test_brent($p1, {a => 0, b => 3}, 2);
test_brent($p1, {a => 3, b => 0}, 2);
test_brent($p1, {a => -1, b => 3, precision => 0.0001}, 2);
test_brent($p1, {a => 0, b => 3, precision => 0.0000001}, 2);
test_brent($p1, {a => 0, b => -100, precision => 0.0000001}, -2);

# what happens if a or b is the root?
test_brent($p1, {a => 2, b => 3}, 2);
is($p1->iterations,0,"a was identified as root at once");
test_brent($p1, {a => 1, b => 2}, 2);
is($p1->iterations,0,"b was identified as root at once");

# a more complicated case
my $p2 = Math::Polynom->new(5 => 5, 3.2 => 4, 0.9 => -2);  # 5*x^5 + 4*x^3.2 - 2*x^0.9
test_brent($p2, {a => 0.5, b => 1, precision => 0.000000000000001}, 0.6161718040343);

eval { test_brent($p2, {a => 0.5, b => 1, precision => 0.000000000000001}, 0.6161718040343); };
ok((!defined $@ || $@ eq ''),"brent() does not fails on polynom 2 with negative guess (newton_raphson would)");

my $p3 = Math::Polynom->new(2 => 1, 1 => -2, 0 => 1); # x^2 -2*x +1
test_brent($p3,{a => 0.5, b => 1},1);
# problem: can't find an interval on which polynom is negative on one side and positive on the other, since always pos
#test_brent($p3,{a => 0, b => -500},1);
#test_brent($p3,{a => 0, b => 99999},1);

# TODO: handle calculation overflows...
my $v;
my $p7 = Math::Polynom->new(5 => 5, 3 => 4, 1 => -2);  # 5*x^5 + 4*x^3 - 2*x
eval { $v = $p7->brent(a => -100000000000000000, b => 999999999999999999999, max_depth => 1); };
ok((defined $@ && $@ =~ /reached maximum number of iterations/),"brent() fails when max_depth reached");
ok($p7->error_message =~ /reached maximum number of iterations/,"\$p7->error_message looks good");
is($p7->error,Math::Polynom::ERROR_MAX_DEPTH,"\$p7->error looks good");
# but we still find the solution if enough depth
test_brent($p7,{a => -100000000000000000, b => 999999999999999999999, max_depth => 150, precision => 0.01,sym => 1}, 0.58893);

# empty polynom error
my $p4 = Math::Polynom->new();
eval { $p4->brent(a => 0, b => 1); };
ok((defined $@ && $@ =~ /empty polynom/),"brent() fails on empty polynom");
ok($p4->error_message =~ /empty polynom/,"\$p4->error_message looks good");
is($p4->error,Math::Polynom::ERROR_EMPTY_POLYNOM,"\$p4->error looks good");



( run in 1.309 second using v1.01-cache-2.11-cpan-71847e10f99 )