File Coverage

blib/lib/Math/Polynomial.pm

Criterion	Covered	Total	%
statement	542	549	98.7
branch	247	256	96.4
condition	102	104	98.0
subroutine	75	76	98.6
pod	56	56	100.0
total	1022	1041	98.1

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.018';
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]" },
    sign_of_coeff => undef,
    plus          => q{ + },
    minus         => q{ - },
    leading_plus  => q{},
    leading_minus => q{- },
    times         => q{ },
    power         => q{^},
    variable      => q{x},
    prefix        => q{(},
    suffix        => q{)},
    wrap          => sub { $_[1] },
);

my $global_string_config = {};

my @tree_defaults = (
    fold_sign     => 0,
    expand_power  => 0,
    group         => sub { $_[0] },
    sign_of_coeff => undef,
    wrap          => sub { $_[1] },
    map {
        my $key = $_;
        ($key => sub { croak "missing parameter: $key" })
    } qw(
        variable
        constant
        negation
        sum
        difference
        product
        power
    ),
);

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

# binary operator wrapper generator
# generates functions to be called via overload:
# - upgrading a non-polynomial operand to a compatible polynomial
# - restoring the original operand order
sub _binary {
    my ($method) = @_;
    return sub {
        my ($this, $that, $reversed) = @_;
        if (!ref($that) || !eval { $that->isa('Math::Polynomial') }) {
            $that = $this->new($that);
        }
        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 number_of_terms {
    my ($this) = @_;
    my $zero   = $this->coeff_zero;
    return scalar grep { $zero != $_ } $this->coeff;
}

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(0 <= $i? $this->coeff($i): ());
    while (0 <= --$i) {
        $result = $result->mul($that)->add_const($this->coeff($i));
    }
    return $result;
}

sub mirror {
    my ($this) = @_;
    my $i = 0;
    return $this->new( map { $i++ & 1? -$_: $_ } $this->coeff );
}

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 is_even {
    my ($this) = @_;
    return $this->is_equal($this->mirror);
}

sub is_odd {
    my ($this) = @_;
    return $this->is_equal($this->mirror->neg);
}

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;
    if ($that->is_monic) {
        undef $hd;
    }
    my @rem  = $this->coeff;
    my $i    = $#rem - @den;
    while (0 <= $i) {
        my $q = pop @rem;
        if (defined $hd) {
            foreach my $r (@rem) {
                $r *= $hd;
            }
        }
        my $j = $i--;
        foreach my $d (@den) {
            $rem[$j++] -= $q * $d;
        }
    }
    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 @binary = ();
    while ($exp > 1) {
        push @binary, 1 & $exp;
        $exp >>= 1;
    }
    my $result = $this;
    while (@binary) {
        $result *= $result;
        $result *= $this if pop @binary;
    }
    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 @binary = ();
    while ($exp > 1) {
        push @binary, 1 & $exp;
        $exp >>= 1;
    }
    my $result = $this;
    while (@binary) {
        $result *= $result;
        $result *= $this if pop @binary;
        $result %= $that;
    }
    return $result;
}

sub exp_mod {
    my ($this, $exp) = @_;
    _check_int($exp);
    return $this->new                   if 0 == $this->degree;
    return $this->new($this->coeff_one) if 0 == $exp;
    my @binary = ();
    while ($exp > 1) {
        push @binary, 1 & $exp;
        $exp >>= 1;
    }
    my $result = $this->new($this->coeff_zero, $this->coeff_one);
    while (@binary) {
        $result *= $result;
        if (pop @binary) {
            local $max_degree;
            $result <<= 1;
        }
        $result %= $this;
    }
    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 inflate {
    my ($this, $exp) = @_;
    my $degree = $this->degree;
    return $this if $degree <= 0;
    _check_int($exp);
    croak 'exponent too large' if
        defined($max_degree) && $degree * $exp > $max_degree;
    return $this->new($this->evaluate($this->coeff_one)) if !$exp;
    my @zeroes = ($this->coeff_zero) x ($exp - 1);
    return $this if !@zeroes;
    my ($const, @coeff) = $this->coeff;
    return $this->new($const, map {@zeroes, $_} @coeff);
}

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) = @_;
    my $n   = $this->coeff_zero;
    my $one = $this->coeff_one;
    return $this->new(
        map { $this->coeff($_) * ($n += $one) } 1..$this->degree
    );
}

sub integrate {
    my ($this, $const) = @_;
    my $zero = $this->coeff_zero;
    my $one  = $this->coeff_one;
    my $n    = $zero;
    if (!defined $const) {
        $const = $zero;
    }
    return $this->new(
        $const,
        map {
            $n += $one;
            my $c = $this->coeff($_);
            $zero == $c? $c: $c / $n
        } 0..$this->degree
    );
}

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

sub _make_ltz {
    my ($config, $zero) = @_;
    return 0 if !$config->{'fold_sign'};
    my $sgn = $config->{'sign_of_coeff'};
    return
        defined($sgn)?
            sub { $sgn->($_[0]) < 0     }:
            sub {        $_[0]  < $zero };
}

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 $zero = $this->coeff_zero;
    my $ltz  = _make_ltz(\%config, $zero);
    my $one  = $this->coeff_one;
    my $with_variable = $config{'with_variable'};
    foreach my $exp ($config{'ascending'}? 0..$max_exp: reverse 0..$max_exp) {
        my $coeff = $this->coeff($exp);

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

        # plus/minus
        if ($ltz && $ltz->($coeff)) {
            $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'},
        $config{'wrap'}->($this, $result),
        $config{'suffix'};
}

sub as_horner_tree {
    my ($this, $params) = @_;
    my %config = (@tree_defaults, %{$params});
    my $exp = $this->degree;
    if ($exp < 0) {
        $exp = 0;
    }
    my $zero = $this->coeff_zero;
    my $ltz  = _make_ltz(\%config, $zero);
    my $coeff = $this->coeff($exp);
    my $first_is_neg = $ltz && $ltz->($coeff);
    if ($first_is_neg) {
        $coeff = -$coeff;
    }
    my $result =
        $exp && $this->coeff_one == $coeff? undef:
        $config{'constant'}->($coeff);
    my $is_sum = 0;
    my $var_is_dynamic = 'CODE' eq ref $config{'variable'};
    my $variable = $var_is_dynamic? undef: $config{'variable'};
    while (0 <= --$exp) {
        $coeff = $this->coeff($exp);
        if ($is_sum) {
            $result = $config{'group'}->($result);
        }
        if ($var_is_dynamic) {
            $variable = $config{'variable'}->();
        }
        my $power = $variable;
        if ($config{'expand_power'}) {
            $is_sum = $zero != $coeff;
        }
        else {
            my $skip = 0;
            while ($zero == $coeff) {
                ++$skip;
                last if $skip > $exp;
                $coeff =  $this->coeff($exp - $skip);
            }
            if ($skip) {
                $exp -= $skip;
                ++$skip if 0 <= $exp;
                $power = $config{'power'}->($variable, $skip) if 1 < $skip;
            }
            $is_sum = 0 <= $exp;
        }
        $result =
            defined($result)? $config{'product'}->($result, $power): $power;
        if ($first_is_neg) {
            $result = $config{'negation'}->($result);
            $first_is_neg = 0;
        }
        if ($is_sum) {
            my $is_neg = $ltz && $ltz->($coeff);
            if ($is_neg) {
                $coeff = -$coeff;
            }
            my $const = $config{'constant'}->($coeff);
            $result = $config{$is_neg? 'difference': 'sum'}->($result, $const);
        }
    }
    if ($first_is_neg) {
        $result = $config{'negation'}->($result);
    }
    return $result;
}

sub as_power_sum_tree {
    my ($this, $params) = @_;
    my %config = (@tree_defaults, %{$params});
    my $exp = $this->degree;
    if ($exp < 0) {
        $exp = 0;
    }
    my $zero = $this->coeff_zero;
    my $ltz  = _make_ltz(\%config, $zero);
    my $one  = $this->coeff_one;
    my $result = undef;
    my $var_is_dynamic = 'CODE' eq ref $config{'variable'};
    my $variable = $var_is_dynamic? undef: $config{'variable'};
    while (0 <= $exp) {
        my $coeff = $this->coeff($exp);

        # skip term?
        next if defined($result) && $zero == $coeff;

        # variable and exponent
        my $term = undef;
        if (0 != $exp) {
            if ($var_is_dynamic) {
                $variable = $config{'variable'}->();
            }
            $term = $variable;
            if (1 != $exp) {
                if ($config{'expand_power'}) {
                    my $todo = $exp-1;
                    for (my $tmp=$exp; $tmp>1; --$tmp) {
                        if ($var_is_dynamic) {
                            $variable = $config{'variable'}->();
                        }
                        $term = $config{'product'}->($term, $variable);
                    }
                }
                else {
                    $term = $config{'power'}->($term, $exp);
                }
            }
        }

        # sign and coefficient
        my $is_neg = $ltz && $ltz->($coeff);
        if ($is_neg) {
            $coeff = -$coeff;
        }
        if (0 == $exp || $one != $coeff) {
            my $const = $config{'constant'}->($coeff);
            $term = 0 == $exp? $const: $config{'product'}->($const, $term);
        }

        # summation
        if (defined $result) {
            $result = $config{$is_neg? 'difference': 'sum'}->($result, $term);
        }
        else {
            $result = $is_neg? $config{'negation'}->($term): $term;
        }
    }
    continue {
        --$exp;
    }
    return $result;
}

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 lcm {
    my $result = shift;
    foreach my $that (@_) {
        my $gcd  = $result->gcd($that);
        $result *= $that->div($gcd) if !$that->is_equal($gcd);
    }
    return $result;
}

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__