# Copyright (c) 2007-2009 by Martin Becker.  All rights reserved.
# This package is free software; you can redistribute it and/or modify it
# under the same terms as Perl itself.
#
# $Id: Polynomial.pm 60 2009-06-11 21:51:47Z demetri $

package Math::Polynomial;

use 5.006;
use strict;
use warnings;
use Carp qw(croak);

require overload;

overload->import(
    q{neg}      => 'neg',
    q{+}        => _binary('add'),
    q{-}        => _binary('sub_'),
    q{*}        => _binary('mul'),
    q{/}        => _binary('div'),
    q{%}        => _binary('mod'),
    q{**}       => _lefty('pow'),
    q{<<}       => _lefty('shift_up'),
    q{>>}       => _lefty('shift_down'),
    q{!}        => 'is_zero',
    q{bool}     => 'is_nonzero',
    q{==}       => _binary('is_equal'),
    q{!=}       => _binary('is_unequal'),
    q{""}       => 'as_string',
    q{fallback} => undef,       # auto-generate trivial substitutions
);

# ----- object definition -----

# Math::Polynomial=ARRAY(...)

# .......... index ..........   # .......... value ..........
use constant F_COEFF  => 0;     # coefficients arrayref, ascending degree
use constant F_ZERO   => 1;     # zero element of coefficient space
use constant F_ONE    => 2;     # unit element of coefficient space
use constant F_CONFIG => 3;     # default stringification configuration
use constant NFIELDS  => 4;

# ----- static data -----

our $VERSION     = '1.002';
our $max_degree  = 10_000;      # limit for power operator

# default values for as_string options
my @string_defaults = (
    ascending     => 0,
    with_variable => 1,
    fold_sign     => 0,
    fold_zero     => 1,
    fold_one      => 1,
    fold_exp_zero => 1,
    fold_exp_one  => 1,
    convert_coeff => sub { "$_[0]" },
    plus          => q{ + },
    minus         => q{ - },
    leading_plus  => q{},
    leading_minus => q{- },
    times         => q{ },
    power         => q{^},
    variable      => q{x},
    prefix        => q{(},
    suffix        => q{)},
);

my $global_string_config = {};

# ----- private/protected subroutines -----

# generic polynomial detection hook (see Math::Polynomial::Generic)
sub _is_generic {
    return 0;
}

# binary operator wrapper generator
# generates functions to be called via overload:
# - upgrading a non-polynomial operand to a compatible polynomial
# - casting a generic operand if appropriate
# - restoring the original operand order
sub _binary {
    my ($method) = @_;
    return sub {
        my ($this, $that, $reversed) = @_;
        if (!ref($that) || !eval { $that->isa('Math::Polynomial') }) {
            if ($this->_is_generic) {
                $that = Math::Polynomial->new($that);
            }
            else {
                $that = $this->new($that);
            }
        }
        if ($this->_is_generic) {
            if (!$that->_is_generic) {
                $this = $this->_cast($that);
            }
        }
        elsif ($that->_is_generic) {
            $that = $that->_cast($this);
        }
        if ($reversed) {
            ($this, $that) = ($that, $this);
        }
        return $this->$method($that);
    };
}

# asymmetrically prototyped binary operator wrapper generator
# generates functions to be called via overload:
# - disallowing reverse order of operands
sub _lefty {
    my ($method) = @_;
    return sub {
        my ($this, $that, $reversed) = @_;
        croak 'wrong operand type' if $reversed;
        return $this->$method($that);
    };
}

# integer argument checker
# - make sure arguments are non-negative integer numbers
sub _check_int {
    foreach my $arg (@_) {
        eval {
            use warnings FATAL => 'all';
            $arg == abs int $arg
        } or croak 'non-negative integer argument expected';
    }
    return;
}

# ----- methods -----

sub new {
    my ($this, @coeff) = @_;
    my $class = ref $this;
    my ($zero, $one, $config);
    if ($class) {
        (undef, $zero, $one, $config) = @{$this};
    }
    else {
        my $sample = @coeff? $coeff[-1]: 1;
        $zero   = $sample - $sample;
        $one    = $sample ** 0;
        $config = undef;
        $class  = $this;
    }
    while (@coeff && $zero == $coeff[-1]) {
        pop @coeff;
    }
    return bless [\@coeff, $zero, $one, $config], $class;
}

sub clone {
    my ($this) = @_;
    return bless [@{$this}], ref $this;
}

sub monomial {
    my ($this, $degree, $coeff) = @_;
    my $zero;
    _check_int($degree);
    croak 'exponent too large'
        if defined($max_degree) && $degree > $max_degree;
    if (ref $this) {
        if (!defined $coeff) {
            $coeff = $this->coeff_one;
        }
        $zero  = $this->coeff_zero;
    }
    else {
        if (!defined $coeff) {
            $coeff = 1;
        }
        $zero  = $coeff - $coeff;
    }
    return $this->new(($zero) x $degree, $coeff);
}

sub from_roots {
    my ($this, @roots) = @_;
    my $one = ref($this)? $this->coeff_one: @roots? $roots[-1] ** 0: 1;
    my $result = $this->new($one);
    foreach my $root (@roots) {
        $result = $result->mul_root($root);
    }
    return $result;
}

sub string_config {
    my ($this, $config) = @_;
    my $have_arg = 2 <= @_;
    if (ref $this) {
        if ($have_arg) {
            $this->[F_CONFIG] = $config;
        }
        else {
            $config = $this->[F_CONFIG];
        }
    }
    else {
        if ($have_arg) {
            # note: do not leave ultimate fallback configuration undefined
            $global_string_config = $config || {};
        }
        else {
            $config = $global_string_config;
        }
    }
    return $config;
}

sub interpolate {
    my ($this, $xvalues, $yvalues) = @_;
    if (
        !ref($xvalues) || !ref($yvalues) || @{$xvalues} != @{$yvalues}
    ) {
        croak 'usage: $q = $p->interpolate([$x1, $x2, ...], [$y1, $y2, ...])';
    }
    return $this->new if !@{$xvalues};
    my @alpha  = @{$yvalues};
    my $result = $this->new($alpha[0]);
    my $aux    = $result->monomial(0);
    my $zero   = $result->coeff_zero;
    for (my $k=1; $k<=$#alpha; ++$k) {
        for (my $j=$#alpha; $j>=$k; --$j) {
            my $dx = $xvalues->[$j] - $xvalues->[$j-$k];
            croak 'x values not disjoint' if $zero == $dx;
            $alpha[$j] = ($alpha[$j] - $alpha[$j-1]) / $dx;
        }
        $aux = $aux->mul_root($xvalues->[$k-1]);
        $result += $aux->mul_const($alpha[$k]);
    }
    return $result;
}

sub coeff {
    my ($this, $degree) = @_;
    if (defined $degree) {
        return
            0 <= $degree && $degree < @{$this->[F_COEFF]}?
                $this->[F_COEFF]->[$degree]:
                $this->[F_ZERO];
    }
    croak 'array context required if called without argument' if !wantarray;
    return @{$this->[F_COEFF]};
}

sub coefficients {
    my ($this) = @_;
    croak 'array context required' if !wantarray;
    return $this->is_zero? ($this->coeff_zero): $this->coeff;
}

sub degree     { return $#{ $_[0]->[F_COEFF] }; }
sub coeff_zero { return $_[0]->[F_ZERO]; }
sub coeff_one  { return $_[0]->[F_ONE]; }

sub proper_degree {
    my ($this) = @_;
    my $degree = $this->degree;
    return 0 <= $degree? $degree: undef;
}

sub evaluate {
    my ($this, $x) = @_;
    my $i = $this->degree;
    my $result = 0 <= $i? $this->coeff($i): $this->coeff_zero;
    while (0 <= --$i) {
        $result = $x * $result + $this->coeff($i);
    }
    return $result;
}

sub nest {
    my ($this, $that) = @_;
    my $i = $this->degree;
    my $result = $that->new;
    while (0 <= $i) {
        $result = $result->mul($that)->add_const($this->coeff($i));
        --$i;
    }
    return $result;
}

sub is_zero {
    my ($this) = @_;
    return $this->degree < 0;
}

sub is_nonzero {
    my ($this) = @_;
    return $this->degree >= 0;
}

sub is_equal {
    my ($this, $that) = @_;
    my $i = $this->degree;
    my $eq = $i == $that->degree;
    while ($eq && 0 <= $i) {
        $eq = $this->coeff($i) == $that->coeff($i);
        --$i;
    }
    return $eq;
}

sub is_unequal {
    my ($this, $that) = @_;
    my $i = $this->degree;
    my $eq = $i == $that->degree;
    while ($eq && 0 <= $i) {
        $eq = $this->coeff($i) == $that->coeff($i);
        --$i;
    }
    return !$eq;
}

sub neg {
    my ($this) = @_;
    return $this if $this->degree < 0;
    return $this->new( map { -$_ } $this->coeff );
}

sub add {
    my ($this, $that) = @_;
    my $this_d = $this->degree;
    my $that_d = $that->degree;
    my $min_d  = $this_d <= $that_d? $this_d: $that_d;
    return $this->new(
        (map { $this->coeff($_) + $that->coeff($_) } 0..$min_d),
        (map { $this->coeff($_) } $that_d+1 .. $this_d),
        (map { $that->coeff($_) } $this_d+1 .. $that_d),
    );
}

sub sub_ {
    my ($this, $that) = @_;
    my $this_d = $this->degree;
    my $that_d = $that->degree;
    my $min_d  = $this_d <= $that_d? $this_d: $that_d;
    return $this->new(
        (map {  $this->coeff($_) - $that->coeff($_) } 0..$min_d),
        (map {  $this->coeff($_) } $that_d+1 .. $this_d),
        (map { -$that->coeff($_) } $this_d+1 .. $that_d),
    );
}

sub mul {
    my ($this, $that) = @_;
    my $this_d = $this->degree;
    return $this if $this_d < 0;
    my $that_d = $that->degree;
    return $this->new(
        map {
            my ($i, $j) = $_ <= $this_d? ($_, 0): ($this_d, $_-$this_d);
            my $sum = $this->coeff($i) * $that->coeff($j);
            while ($i > 0 && $j < $that_d) {
                $sum += $this->coeff(--$i) * $that->coeff(++$j);
            }
            $sum
        } $that_d < 0? (): (0 .. $this_d+$that_d)
    );
}

sub _divmod {
    my ($this, $that) = @_;
    my @den  = $that->coeff;
    @den or croak 'division by zero polynomial';
    my $hd   = pop @den;
    if ($that->is_monic) {
        undef $hd;
    }
    my @rem  = $this->coeff;
    my @quot = ();
    my $i    = $#rem - @den;
    while (0 <= $i) {
        my $q = pop(@rem);
        if (defined $hd) {
            $q /= $hd;
        }
        $quot[$i] = $q;
        my $j = $i--;
        foreach my $d (@den) {
            $rem[$j++] -= $q * $d;
        }
    }
    return (\@quot, \@rem);
}

sub divmod {
    my ($this, $that) = @_;
    croak 'array context required' if !wantarray;
    my ($quot, $rem) = _divmod($this, $that);
    return ($this->new(@{$quot}), @{$quot}? $this->new(@{$rem}): $this);
}

sub div {
    my ($this, $that) = @_;
    my ($quot) = _divmod($this, $that);
    return $this->new(@{$quot});
}

sub mod {
    my ($this, $that) = @_;
    my ($quot, $rem) = _divmod($this, $that);
    return @{$quot}? $this->new(@{$rem}): $this;
}

sub mmod {
    my ($this, $that) = @_;
    my @den  = $that->coeff;
    @den or croak 'division by zero polynomial';
    my $hd   = pop @den;
    my @rem  = $this->coeff;
    my $i    = $#rem - @den;
    while (0 <= $i) {
        my $q = pop @rem;
        my $j = 0;
        while ($j < $i) {
            $rem[$j++] *= $hd;
        }
        foreach my $d (@den) {
            ($rem[$j++] *= $hd) -= $q * $d;
        }
        --$i;
    }
    return $this->new(@rem);
}

sub add_const {
    my ($this, $const) = @_;
    return $this->new($const) if $this->is_zero;
    my $i = 0;
    return $this->new( map { $i++? $_: $_ + $const } $this->coeff );
}

sub sub_const {
    my ($this, $const) = @_;
    return $this->new(-$const) if $this->is_zero;
    my $i = 0;
    return $this->new( map { $i++? $_: $_ - $const } $this->coeff );
}

sub mul_const {
    my ($this, $factor) = @_;
    return $this->new( map { $_ * $factor } $this->coeff );
}

sub div_const {
    my ($this, $divisor) = @_;
    croak 'division by zero' if $this->coeff_zero == $divisor;
    return $this->new( map { $_ / $divisor } $this->coeff );
}

sub mul_root {
    my ($this, $root) = @_;
    return $this->shift_up(1)->sub_($this->mul_const($root));
}

sub _divmod_root {
    my ($this, $root) = @_;
    my $i   = $this->degree;
    my $rem = $this->coeff($i < 0? 0: $i);
    my @quot;
    while (0 <= --$i) {
        $quot[$i] = $rem;
        $rem = $root * $rem + $this->coeff($i);
    }
    return (\@quot, $rem);
}

sub divmod_root {
    my ($this, $root) = @_;
    croak 'array context required' if !wantarray;
    my ($quot, $rem) = _divmod_root($this, $root);
    return ($this->new(@{$quot}), $this->new($rem));
}

sub div_root {
    my ($this, $root, $check) = @_;
    my ($quot, $rem) = _divmod_root($this, $root);
    if ($check && $this->coeff_zero != $rem) {
        croak 'non-zero remainder';
    }
    return $this->new(@{$quot});
}

sub is_monic {
    my ($this) = @_;
    my $degree = $this->degree;
    return 0 <= $degree && $this->coeff_one == $this->coeff($degree);
}

sub monize {
    my ($this) = @_;
    return $this if $this->is_zero || $this->is_monic;
    return $this->div_const($this->coeff($this->degree));
}

sub pow {
    my ($this, $exp) = @_;
    _check_int($exp);
    my $degree = $this->degree;
    return $this->new($this->coeff_one)        if 0 == $exp;
    return $this                                if 0  > $degree;
    return $this->new($this->coeff(0) ** $exp) if 0 == $degree;
    croak 'exponent too large'
        if defined($max_degree) && $degree * $exp > $max_degree;
    my $result = undef;
    while ($exp) {
        if (1 & $exp) {
            $result = defined($result)? $this->mul($result): $this;
        }
        $exp >>= 1 and $this = $this->mul($this);
    }
    return $result;
}

sub pow_mod {
    my ($this, $exp, $that) = @_;
    _check_int($exp);
    $this = $this->mod($that);
    my $this_d = $this->degree;
    return $this->new                          if 0 == $that->degree;
    return $this->new($this->coeff_one)        if 0 == $exp;
    return $this                               if 0  > $this_d;
    return $this->new($this->coeff(0) ** $exp) if 0 == $this_d;
    my $result = undef;
    while ($exp) {
        if (1 & $exp) {
            $result =
                defined($result)? $this->mul($result)->mod($that): $this;
        }
        $exp >>= 1 and $this = $this->mul($this)->mod($that);
    }
    return $result;
}

sub shift_up {
    my ($this, $exp) = @_;
    _check_int($exp);
    croak 'exponent too large' if
        defined($max_degree) && $this->degree + $exp > $max_degree;
    return $this if !$exp;
    return $this->new(($this->coeff_zero) x $exp, $this->coeff);
}

sub shift_down {
    my ($this, $exp) = @_;
    _check_int($exp);
    return $this if !$exp;
    return $this->new( map { $this->coeff($_) } $exp .. $this->degree );
}

sub slice {
    my ($this, $start, $count) = @_;
    _check_int($start, $count);
    my $degree = $this->degree;
    my $end = $start+$count-1;
    if ($degree <= $end) {
        return $this if 0 == $start;
        $end = $degree;
    }
    return $this->new( map { $this->coeff($_) } $start .. $end );
}

sub differentiate {
    my ($this) = @_;
    return $this->new( map { $this->coeff($_) * $_ } 1..$this->degree );
}

sub integrate {
    my ($this, $const) = @_;
    if (!defined $const) {
        $const = $this->coeff_zero;
    }
    return $this->new(
        $const, map { $this->coeff($_) / ($_+1) } 0..$this->degree
    );
}

sub definite_integral {
    my ($this, $lower, $upper) = @_;
    my $ad = $this->integrate;
    return $ad->evaluate($upper) - $ad->evaluate($lower);
}

sub as_string {
    my ($this, $params) = @_;
    my %config = (
        @string_defaults,
        %{$params || $this->string_config || (ref $this)->string_config},
    );
    my $max_exp = $this->degree;
    if ($max_exp < 0) {
        $max_exp = 0;
    }
    my $result = q{};
    my $one = $this->coeff_one;
    foreach my $exp ($config{'ascending'}? 0..$max_exp: reverse 0..$max_exp) {
        my $coeff = $this->coeff($exp);
        my $with_variable = $config{'with_variable'};

        # skip term?
        if (
            $with_variable &&
            $exp < $max_exp &&
            $config{'fold_zero'} &&
            $coeff == $this->coeff_zero
        ) {
            next;
        }

        # plus/minus
        if ($config{'fold_sign'} && $coeff < $this->coeff_zero) {
            $coeff = -$coeff;
            $result .= $config{q[] eq $result? 'leading_minus': 'minus'};
        }
        else{
            $result .= $config{q[] eq $result? 'leading_plus': 'plus'};
        }

        # coefficient
        if (
            !$with_variable ||
            !$config{'fold_one'} ||
            0 == $exp && $config{'fold_exp_zero'} ||
            $one != $coeff
        ) {
            $result .= $config{'convert_coeff'}->($coeff);
            next if !$with_variable;
            if (0 != $exp || !$config{'fold_exp_zero'}) {
                $result .= $config{'times'};
            }
        }

        # variable and exponent
        if (0 != $exp || !$config{'fold_exp_zero'}) {
            $result .= $config{'variable'};
            if (1 != $exp || !$config{'fold_exp_one'}) {
                $result .= $config{'power'} . $exp;
            }
        }
    }
    return join q{}, $config{'prefix'}, $result, $config{'suffix'};
}

sub gcd {
    my ($this, $that, $mod) = @_;
    defined $mod or $mod = 'mod';
    my $mod_op = $this->can($mod);
    $mod_op or croak "no such method: $mod";
    my ($this_d, $that_d) = ($this->degree, $that->degree);
    if ($this_d < $that_d) {
        ($this, $that) = ($that, $this);
        ($this_d, $that_d) = ($that_d, $this_d);
    }
    while (0 <= $that_d) {
        ($this, $that) = ($that, $this->$mod_op($that));
        ($this_d, $that_d) = ($that_d, $that->degree);
        $this_d > $that_d or croak 'bad modulo operator';
    }
    return $this;
}

sub xgcd {
    my ($this, $that) = @_;
    croak 'array context required' if !wantarray;
    my ($d1, $d2) = ($this->new($this->coeff_one), $this->new);
    if ($this->degree < $that->degree) {
        ($this, $that) = ($that, $this);
        ($d1, $d2) = ($d2, $d1);
    }
    my ($m1, $m2) = ($d2, $d1);
    while (!$that->is_zero) {
        my ($div, $mod) = $this->divmod($that);
        ($this, $that) = ($that, $mod);
        ($d1, $d2, $m1, $m2) =
            ($m1, $m2, $d1->sub_($m1->mul($div)), $d2->sub_($m2->mul($div)));
    }
    return ($this, $d1, $d2, $m1, $m2);
}

sub inv_mod {
    my ($this, $that) = @_;
    my ($d, $d2) = ($that->xgcd($this))[0, 2];
    croak 'division by zero polynomial' if $that->is_zero || $d2->is_zero;
    return $d2->div_const($d->coeff($d->degree));
}

1;
__END__

=head1 NAME

Math::Polynomial - Perl class for polynomials in one variable

=head1 VERSION

This documentation refers to version 1.002 of Math::Polynomial.

=head1 SYNOPSIS

  use Math::Polynomial 1.000;

  $p = Math::Polynomial->new(0, -2, 0, 1);    # x^3 - 2 x

  print "p = $p\n";                           # p = (x^3 + -2 x)

  $p->string_config({ fold_sign => 1 });
  print "p = $p\n";                           # p = (x^3 - 2 x)

  $q = $p->new(0, 3, 0, -4, 0, 1);            # x^5 - 4 x^3 + 3 x

  $r = $p ** 2 - $p * $q;                     # arithmetic expression
  $bool = $p == $q;                           # boolean expression

  ($s, $t) = $r->divmod($q);                  # q * s + t = r

  $u = $r->gcd($q);                        # greatest common divisor,
                                           # here: u = 3 x
  $v = $u->monize;                         # v = x
  
  $y = $p->evaluate(0.5);                     # y = p(0.5) = -0.875
  $d = $q->degree;                            # d = degree(q) = 5

  $w = $p->interpolate([0..2], [-1, 0, 3]);   # w(0) = -1, w(1) = 0,
                                              # w(2) = 3

  use Math::Complex;
  $p = Math::Polynomial->new(i, 1+i);         # p(x) = (1+i)*x + i

=head1 DESCRIPTION

Math::Polynomial objects represent polynomials in one variable,
i.e. expressions built with finitely many additions, subtractions
and multiplications of the variable and some constants.  A standard
way of writing down a polynomial in one variable is as a sum of
products of some constant and a power of x, ordered by powers of
x.  The constants in those terms are called coefficients.

The polynomial I<p(x) = 0> is called the zero polynomial.  For
polynomials other than the zero polynomial, the exponent of the
highest power of x with a nonzero coefficient is called the degree
of the polynomial.

New Math::Polynomial objects can be created using a variety of
constructors, or as results of expressions composed from existing objects.
Math::Polynomial objects are immutable with respect to mathematical
properties; all operations on polynomials create and return new objects
rather than modifying anything.

The module works with various types of coefficients, like ordinary
floating point numbers, complex numbers, arbitrary precision
rationals, matrices, elements of finite fields, and lots of others.
All that is required is that the coefficients are either Perl numbers
or objects with suitably overloaded arithmetic operators.  Operations
on polynomials are carried out by reducing them to basic operations
in the domain of their coefficients.

Math::Polynomial objects are implicitly bound to their coefficient
space, which will be inherited when new polynomials are derived
from existing ones, or determined from actual coefficients when
polynomials are created from scratch.  It is the responsibility of
the application not to mix coefficients that cannot be added to or
multiplied by each other.

Note that ordinary Perl numbers used as coefficients have the disadvantage
that rounding errors may lead to undesired effects, such as unexpectedly
non-zero division remainders or failing equality checks.

=head1 CLASS VARIABLES

=over 4

=item I<$VERSION>

C<$VERSION> contains the current version number of the module.  Its
most typical use is the statement:

  use Math::Polynomial 1.000;

This will make sure the module version used is at least 1.000, which
is recommended because previous versions had a different API.

=item I<$max_degree>

C<$max_degree> limits the arguments for the I<pow> and I<shift_up>
operators and the monomial constructor, see L</pow>.  Its default
value is ten thousand.  It can be undefined to disable any size checks.

=back

=head1 CLASS METHODS

=head2 Constructors

=over 4

=item I<new($coeff0, $coeff1, $coeff2, ...)>

C<Math::Polynomial-E<gt>new(@coeff)> creates a new polynomial with
the given coefficients for x to the power of zero, one, two, etc.
For example, C<Math::Polynomial-E<gt>new(7, 1)> creates an object
representing I<p(x) = 7+x>.

Note that coefficients are specified in ascending order of powers
of x.  The degree of the polynomial will be at most I<n-1> if I<n>
coefficients are given, less if one or more of the highest-order
coefficients are zero.

Specifying at least one coefficient (which may be zero) ensures
that the created polynomials use the desired coefficient space.
Without any parameters, I<new> creates a zero polynomial on Perl
numeric values.

=item I<monomial($degree)>

=item I<monomial($degree, $coeff)>

C<Math::Polynomial-E<gt>monomial($degree, $coeff)> creates a
polynomial with $coeff as the coefficient for x to the power of
$degree, and all other coefficients zero.  The degree must be a
non-negative integer number.  If $coeff is omitted, it defaults to
the Perl scalar B<1>.

To prevent accidential excessive memory consumption, C<$degree>
must be at most C<$Math::Polynomial::max_degree>.

=item I<interpolate([$x1, $x2, ...], [$y1, $y2, ...])>

C<Math::Polynomial-E<gt>interpolate(\@x_values, \@y_values)> creates
a Newton/Lagrange interpolation polynomial passing through support
points with the given x- and y-coordinates.  The x-values must be
mutually distinct.  The number of y-values must be equal to the
number of x-values.  For I<n> support points this takes I<O(n**2)>
additions, subtractions, multiplications, divisions and comparisons
each in the coefficient space.  The result will be a polynomial of
degree I<n-1> at most.

Note that with increasing numbers of support points, interpolation
tends to get inaccurate if carried out with limited numerical
precision.  Furthermore, high-degree interpolation polynomials can
oscillate wildly in the neighbourhood of the support points, let
alone elswhere, unless the support point x-values are carefully
chosen.  This is due to the nature of these functions and not a
fault of the module.  A script demonstrating this phenomenon can
be found in the Math-Polynomial examples directory.

=item I<from_roots($x1, $x2, ...)>

C<Math::Polynomial-E<gt>from_roots(@x_values)> creates the monic
polynomial with the given set of roots, i.e. for x-values
I<x1, x2, ..., xn> the product I<(x-x1)*(x-x2)*...*(x-xn)>.  This
is a Polynomial of degree I<n> if I<n> roots are given.  The x-values
can be given in any order and do not need to be distinct.  For I<n>
roots this takes I<O(n*n)> multiplications and subtractions in the
coefficient space.

=back

=head2 Other class methods

Some properties governing default behaviour can be accessed through the
class method I<string_config>.  See L</String Representation>.

=head1 OBJECT METHODS

=head2 Constructors

Each class-level constructor can be used as an object method, too,
i.e. be invoked from an object rather than a class name.  This way,
coefficient space properties are passed on from the invocant object
to the new object, saving some initial coefficient analysis.  Other
properties like per-object stringification settings (explained
below) are inherited likewise.

=over 4

=item I<new($coeff0, $coeff1, $coeff2, ...)>

If C<$p> refers to a Math::Polynomial object, the object method
C<$p-E<gt>new(@coeff)> creates and returns a new polynomial sharing
inheritable properties with C<$p>, but with its own list of coefficients
as specified in C<@coeff>.

=item I<monomial($degree)>

=item I<monomial($degree, $coeff)>

C<$p-E<gt>monomial($degree, $coeff)> creates a monomial like
C<Math::Polynomial-E<gt>monomial($degree, $coeff)>, but sharing
inheritable properties with C<$p>.  If C<$coeff> is omitted it
defaults to the multiplicative unit element of the coefficient space
of C<$p>.

To prevent accidential excessive memory consumption, C<$degree>
must be at most C<$Math::Polynomial::max_degree>.

=item I<interpolate([$x1, $x2, ...], [$y1, $y2, ...])>

C<$p-E<gt>interpolate(\@x_values, \@y_values)> is the object method
usage of the Newton/Lagrange interpolation polynomial constructor
(see above).  It creates a polynomial passing through support points
with the given x- and y-coordinates.  The x-values must be mutually
distinct.  The number of y-values must be equal to the number of
x-values.  All values must belong to the coefficient space of C<$p>.

=item I<from_roots($x1, $x2, ...)>

If C<$p> refers to a Math::Polynomial object, the object method
C<$p-E<gt>from_roots(@x_values)> creates the monic polynomial with
the given set of roots, i.e. for x-values I<x1, x2, ..., xn> the
product I<(x-x1)*(x-x2)*...*(x-xn)>, sharing inheritable properties
with C<$p>.  For I<n> roots this takes I<O(n*n)> multiplications
and subtractions in the coefficient space.

=item I<clone>

C<$p-E<gt>clone> returns a new object equal to C<$p>.  This method
will rarely be needed as Math::Polynomial objects are immutable
except for their stringification configuration (cf. L</string_config>).

=back

=head2 Property Accessors

=over 4

=item I<coefficients>

C<$p-E<gt>coefficients> returns the coefficients of C<$p> in ascending
order of exponents, including zeroes, up to the highest-order
non-zero coefficient.  The result will be a list of I<n+1> coefficients
for polynomials of degree I<n>, or a single zero coefficient for
zero polynomials.

=item I<coeff>

C<$p-E<gt>coeff> returns the coefficients of C<$p> much like
C<$p-E<gt>coefficients>, but for zero polynomials the result will
be an empty list.

(Mnemonic for I<coeff> versus I<coefficients>: Shorter name, shorter
list.)

=item I<coeff($exp)>

C<$p-E<gt>coeff($exp)> returns the coefficient of degree C<$exp>
of C<$p>.  If C<$exp> is less than zero or larger than the degree
of C<$p>, the zero element of the coefficient space is returned.

=item I<coeff_zero>

C<$p-E<gt>coeff_zero> returns the zero element of the coefficient
space of C<$p>, i.e. the neutral element with respect to addition.

=item I<coeff_one>

C<$p-E<gt>coeff_one> returns the multiplicative unit element of the
coefficient space of C<$p>, i.e. the neutral element with respect
to multiplication.

=item I<degree>

C<$p-E<gt>degree> returns B<-1> if C<$p> is a zero polynomial,
otherwise the degree of C<$p>.

=item I<proper_degree>

C<$p-E<gt>proper_degree> returns B<undef> if C<$p> is a zero
polynomial, otherwise the degree of C<$p>.  This can be useful in
order to catch incorrect numerical uses of degrees where zero
polynomials might be involved.

=item I<is_monic>

C<$p-E<gt>is_monic> returns a boolean value which is true if C<$p> is
monic, which means it is not the zero polynomial and its highest-degree
coefficient is equal to one.  Cf. L</monize>.

=back

=head2 Evaluation

=over 4

=item I<evaluate>

C<$p-E<gt>evaluate($x)> computes the value of the polynomial function
given by C<$p> at the position C<$x>.  For polynomials of degree
I<n>, this takes I<n> multiplications and I<n> additions in the
coefficient space.

=back

=head2 Comparison Operators

All comparison operators return boolean results.

Note that there are no order-checking comparisons (E<lt>, E<lt>=, E<gt>,
E<lt>=E<gt>, ...) as neither polynomial nor coefficient spaces in general
need to be ordered spaces.

=over 4

=item C<!>

=item I<is_zero>

C<$p-E<gt>is_zero> or short C<!$p> or C<not $p> checks whether C<$p>
is a zero polynomial.

=item I<is_nonzero>

C<$p-E<gt>is_nonzero> checks whether C<$p> is not a zero polynomial.
This method may be implicitly called if an object is used in boolean
context such as the condition of a while loop or in an expression
with boolean operators such as C<&&> or C<||>.

=item C<==>

=item I<is_equal>

C<$p-E<gt>is_equal($q)> or short C<$p == $q> checks whether C<$p> is
equivalent to C<$q>.  The result is true if both polynomials have the
same degree and the same coefficients.  For polynomials of equal degree
I<n>, this takes at most I<n+1> equality checks in the coefficient space.

Note that I<p == q> implies that I<p(x) == q(x)> for every I<x>,
but the converse implication is not true for some coefficient spaces.

=item C<!=>

=item I<is_unequal>

C<$p-E<gt>is_unequal($q)> or short C<$p != $q> checks whether C<$p>
is not equivalent to C<$q>.  The result is true if the polynomials
have different degree or at least one pair of coefficients of same
degree is different.  For polynomials of equal degree I<n>, this
takes at most I<n+1> equality checks in the coefficient space.

=back

=head2 Arithmetic Operators

=over 4

=item unary C<->

=item I<neg>

C<$p-E<gt>neg> or short C<-$p> calculates the negative of a polynomial.
For a polynomial of degree I<n>, this takes I<n+1> negations in the
coefficient space.

=item C<+>

=item I<add>

C<$p-E<gt>add($q)> or short C<$p + $q> calculates the sum of two polynomials.
For polynomials of degree I<m> and I<n>, this takes I<1+min(m, n)> additions
in the coefficient space.

=item C<->

=item I<sub_>

C<$p-E<gt>sub_($q)> or short C<$p - $q> calculates the difference
of two polynomials.  For polynomials of degree I<m> and I<n>, this
takes I<1+min(m, n)> subtractions in the coefficient space, plus
I<n-m> negations if I<n> is greater than I<m>.

The trailing underscore in the method name may look a bit odd but
will prevent primitive syntax-aware tools from stumbling over
"misplaced" I<sub> keywords.

=item C< *>

=item I<mul>

C<$p-E<gt>mul($q)> or short C<$p * $q> calculates the product of two
polynomials.  For polynomials of degree I<m> and I<n>, this takes
I<(m+1)*(n+1)> multiplications and I<m*n> additions in the coefficient
space.

=item I<divmod>

C<($q, $r) = $p1-E<gt>divmod($p2)> divides a polynomial by another
polynomial and returns the polynomial part of the quotient, and the
remainder.  The second polynomial must not be a zero polynomial.
The remainder is a polynomial of lesser degree than the second
polynomial and satisfies the equation I<$p1 == $p2*$q + $r>.  For
polynomials of degree I<m> and I<n>, I<mE<gt>=n>, this takes I<m+1-n>
divisions, I<(m+1-n)*n> multiplications and I<(m+1-n)*n> subtractions
in the coefficient space.

=item C</>

=item I<div>

C<$p-E<gt>div($q)> or short C<$p / $q> calculates the polynomial part
of the quotient of two polynomials.  This takes the same operations
in the coefficient space as I<divmod>.

=item C<%>

=item I<mod>

C<$p-E<gt>mod($q)> or short C<$p % $q> calculates the remainder from
dividing one polynomial by another.  This takes the same operations
in the coefficient space as I<divmod>.

=item I<mmod>

C<$p-E<gt>mmod($q)> (modified mod) calculates the remainder from
dividing one polynomial by another, multiplied by some constant.
The constant is I<a**d> where I<a> is the highest coefficient of
I<q> and I<d = degree(p) - degree(q) + 1>, if I<degree(p) E<gt> degree(q)>,
otherwise I<d = 0>.  This operation is suitable to substitute I<mod>
in the Euclidean algorithm and can be calculated without division
in the coefficient space.  For polynomials of degree I<m> and I<n>,
I<mE<gt>=n>, this takes I<(m+1-n)*(m+3*n)/2> multiplications and
I<(m+1-n)*n> subtractions in the coefficient space.

=item C< **>

=item I<pow>

C<$p-E<gt>pow($n)> or short C<$p ** $n> calculates a power of
a polynomial.  The exponent C<$n> must be a non-negative integer.
To prevent accidential excessive time and memory consumption, the
degree of the result must be at most C<$Math::Polynomial::max_degree>.
The degree limit can be configured, see L</$max_degree>.  Calculating
the I<n>-th power of a polynomial of degree I<m> takes I<O(m*m*n*n)>
multiplications and additions in the coefficient space.

=item I<pow_mod>

C<$p1-E<gt>pow_mod($n, $p2)> is equivalent to C<($p1 ** $n) % $p2>,
except that the modulo operation is repeatedly applied to intermediate
results in order to keep their degrees small.  The exponent C<$n>
must be a non-negative integer.

=item I<add_const>

=item I<sub_const>

=item I<mul_const>

=item I<div_const>

The arithmetic operations C<add_const>, C<sub_const>, C<mul_const>
and C<div_const> can be used to efficiently add a constant to or
subtract a constant from a polynomial, or multiply or divide a
polynomial by a constant, respectively.

Overloaded arithmetic operators (I<+>, I<->, I<*>, ...) work
with constants in place of polynomial operands, too, by converting
non-polynomial arguments into constant polynomials first.  However,
this usage is both less efficient and less obvious, and therefore
not recommended.

Note that there is no use for a C<mod_const> method, as polynomial
division by a constant always yields a zero remainder.

=item C<E<lt>E<lt>>

=item I<shift_up>

C<$p-E<gt>shift_up($n)> or short C<$p E<lt>E<lt> $n> calculates the
product of a polynomial and a power of x.  The exponent C<$n> must
be a non-negative integer.  To prevent accidential excessive memory
consumption, the degree of the result must be at most
C<$Math::Polynomial::max_degree>.

=item C<E<gt>E<gt>>

=item I<shift_down>

C<$p-E<gt>shift_down($n)> or short C<$p E<gt>E<gt> $n> divides a
polynomial by a power of x and returns the polynomial part of the
result, i.e. discarding negative powers of x.  The exponent C<$n>
must be a non-negative integer.

Shifting up or down is more efficient than multiplication or division
as it does not take any operations in the coefficient space.

=item I<slice>

C<$p-E<gt>slice($m, $n)> is equivalent to:

  $xm = $p->make_monomial($m);
  $xn = $p->make_monomial($n);
  ($p / $xm) % $xn

I.e., it returns a polynomial built from a slice of the coefficients
of the original polynomial starting with degree C<$m>, and at most
C<$n> coefficients.   However, it is more efficient than division
and modulo as it does not perform any operations in the coefficient
space.  The indexes C<$m> and C<$n> must be non-negative integers.

=item I<mul_root>

C<$p-E<gt>mul_root($c)> calculates the product of a polynomial I<p>
and the linear term I<(x - c)>.  The result is a polynomial that
evaluates to zero at the given root I<c>.  For polynomials of degree
I<n>, this takes I<n+1> multiplications and I<n+1> subtractions.

=item I<div_root>

C<$p-E<gt>div_root($c, $check)> with a true value C<$check> checks
whether a polynomial I<p> has a linear factor I<(x - c)> and divides
it by that factor.  If I<c> is not actually a root of the polynomial,
i.e. the division yields a non-zero remainder, an error is reported.
For polynomials of degree I<n>, this takes I<n> multiplications and I<n>
additions.

With a false value for C<$check> or without a second argument, as in
C<$p-E<gt>div_root($c)>, divisibility will not be checked.  The remainder
will simply be discarded.  This may be useful for factoring out a known
root while working under limited precision.

=item I<divmod_root>

C<($q, $r) = $p-E<gt>divmod_root($c)> divides a polynomial I<p>
by a linear factor I<(x - c)> and returns the result and remainder.
The remainder will be a zero polynomial if I<c> is a root of I<p>,
otherwise a constant polynomial.  For polynomials of degree I<n>, this
takes I<n> multiplications and I<n> additions.

This method must be called in array context.

Note that there is no need for a I<mod_root> method, as that would
be indistinguishable from I<evaluate>.

=item I<monize>

C<$p-E<gt>monize> converts an arbitrary non-zero polynomial to a
monic polynomial via division by its highest-order coefficient.
The result will be monic, i.e. with a highest-order coefficient of
one, if the invocant was not the zero polynomial, otherwise the
zero polynomial.  Monization of a non-zero polynomial of degree I<n>
takes I<n> divisions in the coefficient space.

=back

=head2 Assignment Operators

=over 4

=item C<=>

=item C<+=>

=item C<-=>

=item C<*=>

=item C</=>

=item C<%=>

=item C<**=>

=item C<E<lt>E<lt>=>

=item C<E<gt>E<gt>=>

=item C<&&=>

=item C<||=>

As Math::Polynomial objects are immutable with respect to their
arithmetic properties, assignment operators like C<=>, C<+=>, C<-=>,
C<*=> etc. will always replace the object that is being assigned
to, rather than mutate it.  Thus polynomial objects can safely be
"modified" using assignment operators, without side effects on other
variables referencing the same objects.

Note that the short-circuit behaviour of C<&&> and C<||>, which
return the last expression evaluated, implies that C<&&=> and C<||=>
conditionally replace a polynomial by a given expression (another
polynomial, say), not just a boolean value.

=back

=head2 Miscellaneous Operators

=over 4

=item I<nest>

C<$p1-E<gt>nest($p2)> calculates the nested polynomial I<p1(p2(x))>
from two polynomials I<p1(x)> and I<p2(x)>.  For polynomials of
degree I<m> and I<n> this takes I<O(m*m*n*n)> multiplications and
additions in the coefficient space.  The result will be a polynomial
of degree I<m*n> if neither of the polynomials is a zero polynomial,
otherwise a constant or zero polynomial.

=item I<gcd>

C<$p1-E<gt>gcd($p2, $mod)> calculates a greatest common divisor of
two polynomials, using the Euclidean algorithm and the modulo
operator as specified by name.  The C<$mod> parameter is optional
and defaults to C<'mod'>.  With polynomials of degree I<m> and I<n>,
I<mE<gt>=n>, and the default modulo operator I<mod>, this takes
at most I<n> polynomial divisions of decreasing degrees or I<O(m+n)>
divisions and I<O(m*n)> multiplications and subtractions in the
coefficient space.  With the I<mmod> operator, this takes I<O(m*n)>
multiplications and subtractions in the coefficient space.

I<mmod> can have advantages over I<mod> in situations where division
in the coefficient space is much more expensive than multiplication.

=item I<xgcd>

C<($d, $d1, $d2, $m1, $m2) = $p1-E<gt>xgcd($p2)> calculates a
greatest common divisor I<d> and four polynomials I<d1>, I<d2>,
I<m1>, I<m2>, such that I<d = p1*d1 + p2*d2>, I<0 = p1*m1 + p2*m2>,
I<degree(m1*d) = degree(p2)>, I<degree(m2*d) = degree(p1)>, using
the extended Euclidean algorithm.  With polynomials of degree I<m>
and I<n>, I<mE<gt>=n>, this takes at most I<n> polynomial divisions
and I<2*n> polynomial multiplications and subtractions of decreasing
degrees, or, in the coefficient space: I<O(m+n)> divisions and
I<O(m*n)> multiplications and subtractions.

=item I<inv_mod>

C<$q = $p1-E<gt>inv_mod($p2)> calculates the multiplicative pseudo-inverse
I<q> of a polynomial I<p1> modulo another polynomial I<p2>.  I<p2>
must not be a zero polynomial, and I<p1> must not be equivalent to zero
modulo I<p2>.

If I<p1> and I<p2> are relatively prime, the result will be a true
modular inverse, i.e. I<q * p1 % p2> will be the same as I<p2 ** 0>.
In any case, I<q> will be chosen such that I<q * p1 % p2> is the
monic greatest common divisor of I<p1> and I<p2>.

I<inv_mod> takes the same operations in the coefficient space as I<xgcd>
plus I<monize>.

=back

=head2 Calculus Operators

Calculus operators as presented here are meaningful mostly on coefficient
spaces of real or complex numbers.  In particular, they require that
coefficients can be multiplied and divided by integer numbers.

=over 4

=item I<differentiate>

If the coefficient space allows multiplication by Perl integers,
C<$p-E<gt>differentiate> calculates the first derivative of a
polynomial.  For a polynomial of degree I<n>, this takes I<n>
multiplications in the coefficient space.

=item I<integrate>

If the coefficient space allows division by Perl integers,
C<$p-E<gt>integrate> calculates an antiderivative of a polynomial.
The coefficient of degree zero of the result will be zero.
C<$p-E<gt>integrate($c)> does the same but adds the constant C<$c>.
For a polynomial of degree I<n>, both forms of integration take
I<n+1> divisions in the coefficient space.

=item I<definite_integral>

If the coefficient space allows division by Perl integers,
C<$p-E<gt>definite_integral($x1, $x2)> calculates the value of the
definite integral from C<$x1> to C<$x2> over the polynomial function
given by C<$p>.  For real numbers I<x1 E<lt> x2>, this can be
interpreted as the signed area bound by the lines I<x=x1>, I<y=0>,
I<x=x2> and the graph of I<p(x)>, where parts below the x-axis are
regarded as negative.

For a polynomial of degree I<n>, this takes I<n+1> divisions, I<2*n>
multiplications, I<2*n> additions and I<1> subtraction in the
coefficient space.  If you need to calculate more than one definite
integral over the same polynomial function, it is more efficient
to store an antiderivative once (see L</integrate>) and evaluate it
at the different interval limits.  The statement...

  $a = $p->definite_integral($x1, $x2);

... is essentially equivalent to:

  $p_int = $p->integrate;
  $a     = $p_int->evaluate($x2) - $p_int->evaluate($x1);

=back

=head2 String Representation

=over 4

=item C<"">

=item I<as_string>

In string context, Math::Polynomial objects will automatically be
converted to a character string, which is the same as the result
of the I<as_string()> method when called without parameter.

An optional configuration hashref controls many layout aspects of the
string representation.  In the absence of an explicit configuration,
a per-object default configuration is used, and in the absence of that,
a per-class default configuration (see L</string_config>).

Each individual configuration setting has a default value as defined
in the next section.

=item I<string_config>

C<$p-E<gt>string_config($hashref)> sets the per-object default
stringification configuration to C<$hashref>.  C<$hashref> may
be B<undef> to remove a previously set configuration.

C<$p-E<gt>string_config> returns the per-object default stringification
configuration as a reference to a hash, if present, otherwise undef.

C<Math::Polynomial-E<gt>string_config($hashref)> sets the
per-class default stringification configuration to C<$hashref>.
C<Math::Polynomial-E<gt>string_config> returns that configuration.
It should always refer to an existing hash, which may be empty.

A per-object configuration will be propagated to any new objects created
from an object.  Thus it is easy to use consistent settings without
having to touch global parameters.

=back

=head2 Stringification Configuration Options

=over 4

=item ascending

True value: order coefficients from lowest to highest degree;
False value (default): from highest to lowest.

=item with_variable

True value (default): display coefficients together with powers of
the variable; false value: display coefficients alone.  False implies
that I<fold_zero>, I<times>, I<power> and I<variable> will have no
effect.

=item fold_sign

True value: contract the addition symbol and the sign of a negative
value to a single subtraction symbol; false value (default): do not
carry out this kind of replacement.  True is only allowed if the
coefficient space defines a "less than" operator.

=item fold_zero

True value (default): suppress terms with coefficients equal to
zero; false value: do not suppress any terms.  Zero polynomials
are represented with a zero constant term in any case.

=item fold_one

True value (default): suppress coefficients equal to one when
multiplied by a variable power; false value: do not suppress factors
of one.  Note that coefficients very close but not quite equal to
one might be stringified to one without being caught by this rule.

=item fold_exp_zero

True value (default): suppress the variable and the zero exponent
in the constant term; false value: display even the constant term
with a variable power.

=item fold_exp_one

True value (default): suppress the exponent in the term of the
variable to the power of one; false value: display even the linear
term with an exponent.

=item convert_coeff

Code reference specifying a function that takes a coefficient value
and returns a string representing that value.  Default is ordinary
stringification.

=item plus

Addition symbol to put between terms.  Default is a plus character
surrounded by blanks.

=item minus

Subtraction symbol replacing a plus symbol and a negative sign, if
applicable (see I<fold_sign>).  Default is a minus character
surrounded by blanks.

=item leading_plus

Sign symbol to put before the first term unless I<fold_sign> is true
and the coefficient is negative.  Default is an empty string.

=item leading_minus

Sign symbol replacing a negative sign at the first term, if
I<fold_sign> is true.  Default is a minus followed by a blank.

=item times

Multiplication symbol to put between the coefficient and the variable
in each term.  Default is a blank.

=item power

Exponentiation symbol to put between the variable and the exponent
in each term.  Default is a caret (C<^>).

=item variable

Symbol representing the variable.  Default is a lower case x.

=item prefix

Prefix to prepend to the entire polynomial.  Default is a left
parenthesis.

=item suffix

Suffix to append to the entire polynomial.  Default is a right
parenthesis.

=back

=head1 EXAMPLES

  use Math::Polynomial 1.000;

  # simple constructors:
  $p  = Math::Polynomial->new(0, -3, 0, 2);    # (2*x**3 - 3*x)
  $zp = Math::Polynomial->new(0);              # the zero polynomial
  $q  = Math::Polynomial->monomial(3, 2);      # (2*x**3)
  $q  = Math::Polynomial->from_roots(5, 6, 7); # (x-5)*(x-6)*(x-7)

  # Lagrange interpolation:
  $q = Math::Polynomial->interpolate([0..3], [0, -1, 10, 45]);
                                               # (2*x**3 - 3*x)

  # constructors used as object methods:
  $q = $p->new(0, -3, 0, 2);                   # (2*x**3 - 3*x)
  $q = $p->new;                                # the zero polynomial
  $q = $p->monomial(3, 2);                     # (2*x**3)
  $q = $p->monomial(3);                        # (x**3)
  $q = $p->interpolate([0..3], [0, -1, 10, 45]); # (2*x**3 - 3*x)
  $q = $p->from_roots(5, 6, 7);                # (x-5)*(x-6)*(x-7)
  $q = $p->clone;                              # q == p

  # properties
  @coeff = $p->coefficients;                   #    (0, -3, 0, 2)
  @coeff = $p->coeff;                          #    (0, -3, 0, 2)
  $coeff = $p->coeff(0);                       #     0
  $coeff = $p->coeff(1);                       #        -3
  $coeff = $p->coeff(2);                       #            0
  $coeff = $p->coeff(3);                       #               2
  $coeff = $p->coeff(4);                       # 0
  @coeff = $zp->coefficients;                  #    (0)
  @coeff = $zp->coeff;                         #    ()
  $const = $p->coeff_zero;                     # 0
  $const = $p->coeff_one;                      # 1
  $n = $p->degree;                             # 3
  $n = $zp->degree;                            # -1
  $n = $p->proper_degree;                      # 3
  $n = $zp->proper_degree;                     # undef

  # evaluation
  $y = $p->evaluate(4);                        # p(4) == 116

  # comparison
  $bool = !$p;         $bool = $p->is_zero;    # zero polynomial?
  if ($p) ...          if ($p->is_nonzero) ... # not zero polynomial?

  $bool = $p == $q;    $bool = $p->is_equal($q);      # equality
  $bool = $p != $q;    $bool = $p->is_unequal($q);    # inequality

  $p = $q && $r;       $p = $q->is_zero? $q: $r;      # choice
  $p = $q || $r;       $p = $q->is_zero? $r: $q;      # choice

  # arithmetic
  $q = -$p;            $q = $p->neg;           # q(x) == -p(x)
  $q = $p1 + $p2;      $q = $p1->add($p2);     # q(x) == p1(x) + p2(x)
  $q = $p1 - $p2;      $q = $p1->sub_($p2);    # q(x) == p1(x) - p2(x)
  $q = $p1 * $p2;      $q = $p1->mul($p2);     # q(x) == p1(x) * p2(x)
  $q = $p1 / $p2;      $q = $p1->div($p2);     # p1 == q * p2 + r,
                                               #   deg r < deg p2
  $r = $p1 % $p2;      $r = $p1->mod($p2);     # p1 == q * p2 + r,
                                               #   deg r < deg p2
  $q = $p + 3;         $q = $p->add_const(3);  # q(x) == p(x) + 3
  $q = $p - 3;         $q = $p->sub_const(3);  # q(x) == p(x) - 3
  $q = $p * 3;         $q = $p->mul_const(3);  # q(x) == p(x) * 3
  $q = $p / 3;         $q = $p->div_const(3);  # q(x) == p(x) / 3

  $q = $p ** 3;        $q = $p->pow(3);        # q(x) == p(x) ** 3
  $q = $p << 3;        $q = $p->shift_up(3);   # q(x) == p(x) * x**3
  $q = $p >> 3;        $q = $p->shift_down(3); # p == q * x**3 + r,
                                               #   deg r < 3

  $r = $p->slice(0, 3);              # p ==  q * x**3 + r, deg r < 3
  $r = $p->slice(2, 3);              # p == (q * x**3 + r) * x**2 + s,
                                     #   deg r < 3, deg s < 2

  ($q, $r) = $p1->divmod($p2);                 # p1 == q * p2 + r,
                                               #   deg r < deg p2

  $q = $p->mul_root(7);                        # q = p * (x-7)

  $q = $p->div_root(7);                        # p = q * (x-7) + c,
                                               #   c is some constant

  $q = $p->div_root(7, 1);                     # p = q * (x-7)
                                               #   if such q exists,
                                               #   otherwise bang!

  ($q, $r) = $p->divmod_root(7);               # p = q * (x-7) + r,
                                               #   deg r < 1

  $r = $p1->mmod($p2);                         # c * p1 == q * p2 + r,
                                               #   deg r < deg p2,
                                               #   c is a constant

  $q = $p1->nest($p2);                         # q(x) == p1(p2(x))

  $bool = $p->is_monic;                        # whether p is monic

  $q = $p->monize;                             # p(x) == q(x) * c,
                                               #   q is monic or zero

  $r = $p1->pow_mod(3, $p2);                   # p1**3 == q * p2 + r,
                                               #   deg r < deg p2

  # greatest common divisor
  $d = $p1->gcd($p2);                          # p1 == q1 * d,
                                               # p2 == q2 * d,
                                               # deg d is maximal

  $d = $p1->gcd($p2, 'mmod');                  # like gcd, but with
                                               # modified modulo op

  # extended Euclidean algorithm
  ($d, $d1, $d2, $m1, $m2) = $p1->xgcd($p2);   # p1 == q1 * d,
                                               # p2 == q2 * d,
                                               # deg d is maximal,
                                               # d == p1*d1 + p2*d2,
                                               # 0 == p1*m1 + p2*m2

  $r = $p1->inv_mod($p2);                      # p1 * r == p2 * q + d,
                                               #   deg r < deg p2,
                                               #   d is monic gcd of
                                               #   p1 and p2

  # calculus
  $q = $p->differentiate;             # q(x) == d/dx p(x)
  $q = $p->integrate;                 # d/dx q(x) == p(x), q(0) == 0
  $q = $p->integrate(10);             # d/dx q(x) == p(x), q(0) == 10
  $a = $p->definite_integral(1, 2);   # Integral from 1 to 2, p(x) dx

  # configurable string representation
  $config = {
    ascending     => 0,
    with_variable => 1,
    fold_sign     => 0,
    fold_zero     => 1,
    fold_one      => 1,
    fold_exp_zero => 1,
    fold_exp_one  => 1,
    convert_coeff => sub { "$_[0]" },
    plus          => q{ + },
    minus         => q{ - },
    leading_plus  => q{},
    leading_minus => q{- },
    times         => q{ },
    power         => q{^},
    variable      => q{x},
    prefix        => q{(},
    suffix        => q{)},
  };
  $str = "$p";                                 # '(2 x^3 + -3 x)'
  $str = $p->as_string();                      # '(2 x^3 + -3 x)'
  $str = $p->as_string($config);               # '(2 x^3 + -3 x)'

  $str = $p->as_string({fold_sign => 1});      # '(2 x^3 - 3 x)'

  $p->string_config({fold_sign => 1});
  $str = "$p";                                 # '(2 x^3 - 3 x)'
  $str = $p->as_string;                        # '(2 x^3 - 3 x)'

  $p->string_config(undef);
  $str = "$p";                                 # '(2 x^3 + -3 x)'

  Math::Polynomial->string_config({fold_sign => 1});
  $str = "$p";                                 # '(2 x^3 - 3 x)'

  $config = $p->string_config;                 # undef
  $config = Math::Polynomial->string_config;   # {fold_sign => 1}

  $config = {
    with_variable => 0,
    ascending     => 1,
    plus          => q{, },
  };
  $str = $p->as_string($config);               # '(0, -3, 0, 2)'

  # other customizations

  $q = do {
    local Math::Polynomial::max_degree;        # override limitation
    $p->monomial(100_000)                      # x^100000
  };

  # examples of other coefficient spaces

  use Math::Complex;
  $c0 = Math::Complex->make(0, 3);
  $c1 = Math::Complex->make(2, 1);
  $p = Math::Polynomial->new($c0, $c1);        # p(x) == (2+i)*x + 3i
  $p->string_config({
    convert_coeff => sub {
      my $coeff = "$_[0]";
      if ($coeff =~ /[+-]/) {
        $coeff = "($coeff)";
      }
      return $coeff;
    },
  });
  print "$p\n";                                # prints ((2+i) x + 3i)

  use Math::BigRat;
  $c0 = Math::BigRat->new('-1/2');
  $c1 = Math::BigRat->new('0');
  $c2 = Math::BigRat->new('3/2');
  $p = Math::Polynomial->new($c0, $c1, $c2);  # p(x) == 3/2*x**2 - 1/2

=head1 DIAGNOSTICS

Currently, Math::Polynomial does not thoroughly check coefficient
values for sanity.  Incompatible coefficients might trigger warnings
or error messages from the coefficient class implementation, blaming
Math::Polynomial to operate on incompatible operands.  Unwisely
chosen coefficient objects, lacking overrides for arithmetic
operations, might even silently be used with their memory address
acting as an integer value.

Some types of wrong usage, however, are diagnosed and will trigger
one of the error messages listed below.  All of them terminate
program execution unless they are trapped in an eval block.

=over 4

=item array context required

A method designed to return more than one value, like I<divmod> or
I<xgcd>, was not called in array context.  The results of these
methods have to be assigned to a list of values rather than a single
value.

=item array context required if called without argument

The I<coeff> method was called without arguments, so as to return
a list of coefficients, but not in array context.  The list of
coefficients should be assigned to an array.  To retrieve a single
coefficient, I<coeff> should be called with a degree argument.

=item bad modulo operator

A modulo operator was passed to I<gcd> and turned out to violate
rules required by such an operator, like constraints on the degree
of returned polynomials.

=item division by zero

The I<div_const> method was called with a zero argument.  Coefficient
spaces are supposed to not allow division by zero, therefore
Math::Polynomial tries to safeguard against it.

=item division by zero polynomial

Some kind of polynomial division, like I<div>, I<mod>, I<divmod>,
I<inv_mod> or an expression with C</> or C<%> was attempted with a
zero polynomial acting as denominator or modulus.  Division by a
zero polynomial is not defined.

=item exponent too large

One of the methods I<pow>, I<shift_up> or I<monomial> was called
with arguments that would lead to a result with an excessively high
degree.  You can tweak the class variable I<$max_degree> to change
the actual limit or disable the check altogether.  Calculations
involving large polynomials can consume a lot of memory and CPU
time.  Exponent sanity checks help to avoid that from happening by
accident.

=item no such method: %s

A modulo operator name was passed to I<gcd> which is not actually
the name of a method.

=item non-negative integer argument expected

One of the methods expecting non-negative integer arguments, like
I<pow>, I<pow_mod>, I<shift_up>, I<shift_down>, I<slice> or
I<monomial>, got something else instead.

=item non-zero remainder

The I<div_root> method was called with a first argument that was not
actually a root of the polynomial and a true value as second argument,
forcing a check for divisibility.

=item usage: %s

A method designed to be called in a certain manner with certain types
of arguments got not what it expected.

For example, I<interpolate> takes two references of arrays of equal
length.

Usage messages give an example of the expected calling syntax.

=item wrong operand type

An arithmetic expression used an operation not defined for polynomial
objects, such as a power with a polynomial exponent.

=item x values not disjoint

The I<interpolate> method was called with at least two equal x-values.
Support points for this kind of interpolation must have distinct x-values.

=back

=head1 EXTENSION GUIDELINES

Math::Polynomial can be extended in different ways.  Subclasses may
restrict coefficient spaces to facilitate certain application domains
such as numerical analysis or channel coding or symbolic algebra.

Special polynomials (such as Legendre, Chebyshev, Gegenbauer, Jacobi,
Hermite, Laguerre polynomials etc.) may be provided by modules
adding not much more than a bunch of funnily named constructors to
the base class.

Other modules may implement algorithms employing polynomials (such
as approximation techniques) or analyzing them (such as root-finding
techniques).

Yet another set of modules may provide alternative implementations
optimized for special cases such as sparse polynomials, or taking
benefit from specialized external math libraries.

Multivariate polynomials, finally, will presumably live in a class
that is neither a base class nor a subclass of Math::Polynomial,
as both subjects do not have enough in common to suggest otherwise.

More detailed information as well as a couple of examples are
supposed to be added in an upcoming release.

=head1 DEPENDENCIES

This version of Math::Polynomial requires these other modules and
libraries to run:

=over 4

=item *

perl version 5.6.0 or higher

=item *

overload (usually bundled with perl)

=item *

Carp (usually bundled with perl)

=back

Additional requirements to run the test suite are:

=over 4

=item *

Test (usually bundled with perl)

=item *

Math::Complex (usually bundled with perl)

=item *

File::Spec (usually bundled with perl)

=item *

File::Basename (usually bundled with perl)

=item *

FindBin (usually bundled with perl)

=back

Recommended modules for increased functionality are:

=over 4

=item *

Math::Complex (usually bundled with perl)

=item *

Math::BigRat

=back

=head1 BUGS AND LIMITATIONS

At the time of release, there were no known unresolved issues with
this module.  Bug reports and suggestions are welcome -- please
submit them through the CPAN RT,
L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Math-Polynomial>.

=head1 SEE ALSO

Part of this distribution:

=over 4

=item *

Example scripts distributed with Math::Polynomial.

Scripts in the examples directory cover module version differences,
interpolation, rational coefficients, pretty-printing, and more.

Scripts in the test suite are less wordily documented, but should have at
least one usage example for each of even the most exotic features.

=item *

I<Math::Polynomial::Generic> -
an experimental interface extension with some syntactical sugar
coating Math::Polynomial.

=back

Modules planned for release, in various states of completion:

=over 4

=item *

I<Math::Polynomial::Sparse> -
an alternative implementation optimized for polynomials with lots
of zero coefficients.

=item *

I<Math::Polynomial::Orthogonal> -
an extension providing constructors for many well-known kinds of
orthogonal polynomials.

=item *

I<Math::Polynomial::Roots::DurandKerner> -
an extension implementing the Weierstrass-(Durand-Kerner) numerical
root-finding algorithm.

=item *

I<Math::Polynomial::Multivariate> -
a complementary module dealing with polynomials in more than one variable.

=item *

I<Math::Polynomial::Parse> -
an extension providing methods to create polynomial objects from strings.

=back

Other Modules:

=over 4

=item *

I<PDL> -
the Perl Data Language.

=item *

Math::GMP -
Interface to the GMP arbitrary precision integer math library.

=item *

Math::Pari -
Interface to the Pari math library.

=back

General Information:

=over 4

=item *

Pages in category I<Polynomials> of Wikipedia.
L<http://en.wikipedia.org/wiki/Category%3APolynomials>

=item *

Weisstein, Eric W. "Polynomial."
From MathWorld--A Wolfram Web Resource.
L<http://mathworld.wolfram.com/Polynomial.html>

=item *

Knuth, Donald E. "The Art of Computer Programming, Volume 2:
Seminumerical Algorithms", 3rd ed., 1998, Addison-Wesley Professional.

=back

=head1 AUTHOR

Martin Becker, E<lt>becker-cpan-mp@cozap.comE<gt>

=head1 ACKNOWLEDGEMENTS

Math::Polynomial was inspired by a module of the same name written and
maintained by Mats Kindahl, E<lt>mats@kindahl.netE<gt>, 1997-2007,
who kindly passed it on to the author for maintenance and improvement.

Additional suggestions and bug reports came from Thorsten Dahlheimer
and Kevin Ryde.

=head1 LICENSE AND COPYRIGHT

Copyright (c) 2007-2009 by Martin Becker.  All rights reserved.

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.6.0 or,
at your option, any later version of Perl 5 you may have available.

=head1 DISCLAIMER OF WARRANTY

This library is distributed in the hope that it will be useful,
but without any warranty; without even the implied warranty of
merchantability or fitness for a particular purpose.

=cut