Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
245 lines (208 sloc) 6.57 KB
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;