``````package Math::GF;
use strict;
use warnings;
{ our \$VERSION = '0.004'; }

use Moo;
use Ouch;
use Math::GF::Zn;

use constant MARGIN => 1.1;

has order          => (is => 'ro');
has p              => (is => 'ro');
has n              => (is => 'ro');
has order_is_prime => (is => 'ro');
has element_class  => (is => 'ro');

# The following are used only for extension fields
has sum_table      => (is => 'ro');
has prod_table     => (is => 'ro');

# neutral element for "+" operation
sub additive_neutral { return \$_[0]->e(0) }

# factory method to create "all" elements in the field
sub all {
my \$self   = shift;
my \$eclass = \$self->element_class;
my \$order  = \$self->order;
map { \$eclass->new(field => \$self, v => \$_) } 0 .. (\$order - 1);
} ## end sub all

# import a handy factory method into caller's package
sub import_builder {
my (\$package, \$order) = splice @_, 0, 2;
my %args = (@_ && ref(\$_[0]) eq 'HASH') ? %{\$_[0]} : @_;

my \$field = \$package->new(order => \$order);
my \$builder = sub { return \$field->e(@_) };
my \$callpkg = caller(\$args{level} // 0);
my \$name = \$args{name} // (
\$field->order_is_prime
? "GF_\$order"
: join('_', 'GF', \$field->p, \$field->n)
);
no strict 'refs';
*{\$callpkg . '::' . \$name} = \$builder;
return;
} ## end sub import_builder

# factory method to create "e"lements of the field
sub e {
my \$self = shift;
my \$ec = \$self->element_class;
return \$ec->new(field => \$self, v => \$_[0]) unless wantarray;
return map { \$ec->new(field => \$self, v => \$_) } @_;
}

# neutral element for "*" operation
sub multiplicative_neutral { return \$_[0]->e(1) }

sub BUILDARGS {
my (\$class, %args) = @_;

ouch 500, 'missing order' unless exists \$args{order};
my \$order = \$args{order};
ouch 500, 'undefined order' unless defined \$order;
ouch 500, 'order MUST be integer and greater than 1'
unless \$order =~ m{\A(?: [2-9] | [1-9]\d+)\z}mxs;

my (\$p, \$n) = __prime_power_decomposition(\$order);
ouch 500, 'order MUST be a power of a prime'
unless defined \$p;
\$args{p}              = \$p;
\$args{n}              = \$n;
\$args{order_is_prime} = (\$n == 1);
if (\$n == 1) {
\$args{order_is_prime} = 1;
\$args{element_class} = 'Math::GF::Zn';
delete @args{qw< sum_table prod_table >};
}
else {
\$args{order_is_prime} = 0;
\$args{element_class} = 'Math::GF::Extension';
@args{qw< sum_table prod_table >} = __tables(\$order);
require Math::GF::Extension;
}

return {%args};
} ## end sub BUILDARGS

sub __tables {
my \$order = shift;

# Get the basic subfield
my (\$p, \$n) = __prime_power_decomposition(\$order);
my \$Zp = Math::GF->new(order => \$p);
my @Zp_all = \$Zp->all;
my (\$zero, \$one) = (\$Zp->additive_neutral, \$Zp->multiplicative_neutral);

my \$pirr = __get_irreducible_polynomial(\$Zp, \$n);
my \$polys = __generate_polynomials(\$Zp, \$n);
my %id_for = map {; "\$polys->[\$_]" => \$_ } 0 .. \$#\$polys;

my (@sum, @prod);
for my \$i (0 .. \$#\$polys) {
my \$I = \$polys->[\$i];
push @sum, \my @ts;
push @prod, \my @tp;
for my \$j (0 .. \$i) {
my \$J = \$polys->[\$j];
my \$sum = (\$I + \$J);
push @ts, \$id_for{"\$sum"};
my \$prod = (\$I * \$J) % \$pirr;
push @tp, \$id_for{"\$prod"};
}
}

return (\@sum, \@prod);
}

sub __generate_polynomials {
my (\$field, \$degree) = @_;
ouch 500, 'irreducible polynomial search only over Zn field'
unless \$field->order_is_prime;
my \$zero = \$field->additive_neutral;
my \$one  = \$field->multiplicative_neutral;

my @coeffs = (\$zero) x (\$degree + 1); # last one as a flag
my @retval;
while (\$coeffs[-1] == \$zero) {
push @retval, Math::Polynomial->new(@coeffs);
for (@coeffs) {
\$_ = \$_ + \$one;
last unless \$_ == \$zero;
}
}
return \@retval;
}

sub __get_irreducible_polynomial {
my (\$field, \$degree) = @_;
ouch 500, 'irreducible polynomial search only over Zn field'
unless \$field->order_is_prime;

my \$zero = \$field->additive_neutral;
my \$one  = \$field->multiplicative_neutral;
require Math::Polynomial;
my @coeffs = (\$one, ((\$zero) x (\$degree - 1)), \$one);
while (\$coeffs[-1] == \$one) {
my \$poly = Math::Polynomial->new(@coeffs);
return \$poly if __rabin_irreducibility_test(\$poly);
for (@coeffs) {
\$_ = \$_ + \$one;
last unless \$_ == \$zero; # wrapped up
}
}
ouch 500, "no monic irreducibile polynomial!"; # never happens
}

sub __to_poly {
my (\$x, \$n) = @_;
my @coeffs;
while (\$x) {
push @coeffs, \$x % \$n;
\$x = (\$x - \$coeffs[-1]) / \$n;
}
push @coeffs, 0 unless @coeffs;
return Z_poly(\$n, @coeffs);
}

sub __rabin_irreducibility_test {
my \$f    = shift;
my \$n    = \$f->degree;
my \$one  = \$f->coeff_one;
my \$pone = Math::Polynomial->monomial(0, \$one);
my \$x    = Math::Polynomial->monomial(1, \$one);
my \$q    = \$one->n;
my \$ps   = __prime_divisors_of(\$n);

for my \$pi (@\$ps) {
my \$ni  = \$n / \$pi;
my \$qni = \$q**\$ni;
my \$h = (Math::Polynomial->monomial(\$qni, \$one) - \$x) % \$f;
my \$g = \$h->gcd(\$f, 'mod');
#return if \$g != \$pone;
return if \$g->degree > 1;
} ## end for my \$pi (@\$ps)
my \$t = (Math::Polynomial->monomial(\$q**\$n, \$one) - \$x) % \$f;
return \$t->degree == -1;
} ## end sub rabin_irreducibility_test

sub __prime_power_decomposition {
my \$x = shift;
return if \$x < 2;
return (\$x, 1) if \$x < 4;

my \$p = __prime_divisors_of(\$x, 'first only please');
return (\$x, 1) if \$x == \$p;    # \$x is prime

my \$e = 0;
while (\$x > 1) {
return if \$x % \$p;         # not the only divisor!
\$x /= \$p;
++\$e;
}
return (\$p, \$e);
} ## end sub __prime_power_decomposition

sub __prime_divisors_of {
my (\$n, \$first_only) = @_;
my @retval;

return if \$n < 2;

for my \$p (2, 3) { # handle cases for 2 and 3 first
next if \$n % \$p;
return \$p if \$first_only;
push @retval, \$p;
\$n /= \$p until \$n % \$p;
}

my \$p = 5; # tentative divisor, will increase through iterations
my \$top = int(sqrt(\$n) + MARGIN); # top attempt for divisor
my \$d = 2; # increase for \$p, alternates between 4 and 2
while (\$p <= \$top) {
if (\$n % \$p == 0) {
return \$p if \$first_only;
unshift @retval, \$p;
\$n /= \$p until \$n % \$p;
\$top = int(sqrt(\$n) + MARGIN);
}
\$p += \$d;
\$d = (\$d == 2) ? 4 : 2;
} ## end while (\$n > 1)

# exited with \$n left as a prime... maybe
return \$n if \$first_only; # always in this case
push @retval, \$n if \$n > 1;

return \@retval;
} ## end sub prime_divisors_of

1;
``````