File Coverage

blib/lib/Math/Polynomial/ModInt.pm

Criterion	Covered	Total	%
statement	199	199	100.0
branch	72	72	100.0
condition	27	27	100.0
subroutine	37	37	100.0
pod	17	17	100.0
total	352	352	100.0

# Copyright (c) 2013-2022 by Martin Becker.  This package is free software,
# licensed under The Artistic License 2.0 (GPL compatible).

package Math::Polynomial::ModInt;

use 5.006;
use strict;
use warnings;
use Math::BigInt try => 'GMP';
use Math::ModInt 0.012;
use Math::ModInt::Event qw(UndefinedResult);
use Math::Polynomial 1.015;
use Carp qw(croak);

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

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

# .............. index ..............   # .......... value ..........

use constant _OFFSET  => Math::Polynomial::_NFIELDS;
use constant _F_INDEX => _OFFSET + 0;    # ordinal number i, 0 <= i
use constant _NFIELDS => _OFFSET + 1;

# ----- class data -----

BEGIN {
    require Exporter;
    our @ISA       = qw(Math::Polynomial Exporter);
    our @EXPORT_OK = qw(modpoly);
    our $VERSION   = '0.005';
    our @CARP_NOT  = qw(Math::ModInt::Event::Trap Math::Polynomial);
}

my $default_string_config = {
    convert_coeff => sub { $_[0]->residue },
    times         => q[*],
    wrap          => \&_wrap,
};

my $lifted_string_config = {
    times         => q[*],
    fold_sign     => 1,
};

my $ipol = Math::Polynomial->new(Math::BigInt->new);
$ipol->string_config($lifted_string_config);

# ----- private subroutines -----

# event catcher
sub _nonprime_modulus {
    croak 'modulus not prime';
}

# event catcher
sub _no_inverse {
    croak 'undefined inverse';
}

# constructor diagnostic
sub _modulus_too_small {
    croak 'modulus must be greater than one';
}

# wrapper for modular inverse to bail out early if not in a field
sub _inverse {
    my ($element) = @_;
    my $trap = UndefinedResult->trap(\&_no_inverse);
    return $element->inverse;
}

# Does a set of coefficient vectors have a zero linear combination?
# We transform and include them one by one into a matrix in echelon form.
# Argument is an iterator so we can short-cut generating superfluous vectors.
# Supplied vectors may be modified.
sub _linearly_dependent {
    my ($it) = @_;
    my $trap = UndefinedResult->trap(\&_nonprime_modulus);
    my @ech = ();
    while (defined(my $vec = $it->())) {
        my $ex = 0;
        for (; $ex < @ech; ++$ex) {
            my $evec = $ech[$ex];
            last if @{$evec} < @{$vec};
            if (@{$evec} == @{$vec}) {
                my $i = $#{$vec};
                return 1 if !$i;
                my $w = pop(@{$vec}) / $evec->[$i];
                while ($i-- > 0) {
                    $vec->[$i] -= $evec->[$i] * $w;
                }
                while (@{$vec} && !$vec->[-1]) {
                    pop(@{$vec});
                }
            }
        }
        return 1 if !@{$vec};
        splice @ech, $ex, 0, $vec;
    }
    return 0;
}

# ----- private methods -----

# wrapper to decorate stringified polynomial
sub _wrap {
    my ($this, $text) = @_;
    my $modulus = $this->modulus;
    return "$text (mod $modulus)";
}

# ----- protected methods -----

# constructor with index argument
sub _xnew {
    my $this  = shift;
    my $index = shift;
    my $poly  = $this->new(@_);
    $poly->[_F_INDEX] = $index;
    return $poly;
}

# ----- overridden public methods -----

sub new {
    my $this = shift;
    if (!@_ && !ref $this) {
        croak 'insufficient arguments';
    }
    if (grep {$_->is_undefined} @_) {
        _nonprime_modulus();
    }
    if (grep {$_->modulus <= 1} @_) {
        _modulus_too_small();
    }
    return $this->SUPER::new(@_);
}

sub string_config {
    my $this = shift;
    return $this->SUPER::string_config(@_) if ref $this;
    ($default_string_config) = @_          if @_;
    return $default_string_config;
}

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

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

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

sub monize {
    my ($this) = @_;
    my $degree = $this->degree;
    return $this if $degree < 0;
    my $leader = $this->coeff($degree);
    return $this if $leader->residue == 1;
    return $this->mul_const(_inverse($leader));
}

# ----- public subroutine -----

sub modpoly { __PACKAGE__->from_index(@_) }

# ----- class-specific public methods -----

sub from_index {
    my ($this, $index, $modulus) = @_;
    my $zero;
    if (defined $modulus) {
        $modulus > 1 || _modulus_too_small();
        $zero = Math::ModInt::mod(0, $modulus);
    }
    elsif (ref $this) {
        $zero = $this->coeff_zero;
        $modulus = $zero->modulus;
    }
    else {
        croak('usage error: modulus parameter missing');
    }
    my @coeff = ();
    my $q = $index;
    while ($q > 0) {
        ($q, my $r) = $zero->new2($q);
        push @coeff, $r;
    }
    return $this->_xnew(@coeff? ($index, @coeff): (0, $zero));
}

sub from_int_poly {
    my ($this, $poly, $modulus) = @_;
    if (!defined $modulus) {
        if (!ref $this) {
            croak('usage error: modulus parameter missing');
        }
        $modulus = $this->modulus;
    }
    my @coeff = map { Math::ModInt::mod($_, $modulus) } $poly->coefficients;
    return $this->new(@coeff);
}

sub modulus {
    my ($this) = @_;
    return $this->coeff_zero->modulus;
}

sub index {
    my ($this) = @_;
    my $base  = $this->modulus;
    my $index = $this->[_F_INDEX];
    if (!defined $index) {
        $index = Math::BigInt->new;
        foreach my $c (reverse $this->coeff) {
            $index->bmul($base)->badd($c->residue);
        }
        $this->[_F_INDEX] = $index;
    }
    return $index;
}

sub number_of_terms { scalar grep { $_->is_not_zero } $_[0]->coeff }

sub lift       { $ipol->new(map { $_->residue          } $_[0]->coefficients) }

sub centerlift { $ipol->new(map { $_->centered_residue } $_[0]->coefficients) }

sub lambda_reduce {
    my ($this, $lambda) = @_;
    my $degree = $this->degree;
    return $this if $degree <= $lambda;
    my @coeff = map { $this->coeff($_) } 0 .. $lambda;
    for (my $i = 1, my $n = $lambda+1; $n <= $degree; ++$i, ++$n) {
        $coeff[$i] += $this->coeff($n);
        $i = 0 if $i == $lambda;
    }
    return $this->new(@coeff);
}

sub first_root {
    my ($this) = @_;
    my $zero = $this->coeff_zero;
    my $lsc  = $this->coeff(0);
    return $zero if $lsc->is_zero;
    if ($this->degree == 1) {
        my $mscr = $this->coeff(1)->centered_residue;
        return -$lsc if $mscr ==  1;
        return  $lsc if $mscr == -1;
    }
    foreach my $n (1 .. $zero->modulus - 1) {
        my $root = $zero->new($n);
        return $root if $this->evaluate($root)->is_zero;
    }
    return undef;
}

# implementation restriction: defined for prime moduli only
sub is_irreducible {
    my ($this) = @_;
    my $n = $this->degree;
    return 0 if $n <= 0;
    return 1 if $n == 1;
    return 0 if !$this->coeff(0);
    return 0 if $this->gcd($this->differentiate)->degree > 0;
    my $p  = $this->modulus;
    # optimization: O(p) zero search only for small p or very large n
    if ($p <= 43 || log($p - 20) <= $n * 0.24 + 2.68) {
        my $rp = 2 < $p && $p <= $n? $this->lambda_reduce($p-1): $this;
        return 0 if defined $rp->first_root;
        return 1 if $n <= 3;
    }
    # Berlekamp rank < $n - 1?
    my $xp = $this->exp_mod($p);
    my $aj = $xp;
    my $bj = $this->monomial(1);
    my $j  = 0;
    return 0 if _linearly_dependent(
        sub {
            return undef if $j >= $n - 1;
            if ($j++) {
                $aj = $aj * $xp % $this;
                $bj <<= 1;
            }
            return [($aj - $bj)->coeff];
        }
    );
    return 1;
}

1;

__END__