Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Rewrite chunks of BigPolynomial

  • Loading branch information...
commit f6a2bc30e14f34f88da4b8f77c47bd3e25d0d260 1 parent 91fceff
@danaj danaj authored
Showing with 127 additions and 107 deletions.
  1. +2 −0  Changes
  2. +12 −12 lib/Math/Primality/AKS.pm
  3. +113 −95 lib/Math/Primality/BigPolynomial.pm
View
2  Changes
@@ -9,6 +9,8 @@
Change AKS to the latest standard implementation (v6 of the paper).
+ Rewrite big parts of BigPolynomial.
+
=head2 0.08
Use a small GCD in the prime_count loop, making it much faster.
View
24 lib/Math/Primality/AKS.pm
@@ -128,22 +128,22 @@ sub is_aks_prime($) {
my $intr = Rmpz_get_ui($r);
debug "Running poly check, r = $intr polylimit = $polylimit\n";
+ my $final_size = Math::GMPz->new(0);
+ Rmpz_mod($final_size, $n, $r);
+ my $compare = Math::Primality::BigPolynomial->new();
+ $compare->setCoef(Math::GMPz->new(1), $final_size);
+
for(my $a = 1; Rmpz_cmp_ui($polylimit, $a) >= 0; $a++) {
debug "Checking at $a / $polylimit\n";
- my $final_size = Math::GMPz->new(0);
- Rmpz_mod($final_size, $n, $r);
- my $compare = Math::Primality::BigPolynomial->new(Rmpz_get_ui($final_size));
- $compare->setCoef(Math::GMPz->new(1), $final_size);
- $compare->setCoef(Math::GMPz->new($a), 0);
- my $res = Math::Primality::BigPolynomial->new($intr);
- my $base = Math::Primality::BigPolynomial->new(1);
- $base->setCoef(Math::GMPz->new(0), $a);
- $base->setCoef(Math::GMPz->new(1), 1);
-
- Math::Primality::BigPolynomial::mpz_poly_mod_power($res, $base, $n, $n, $intr);
+ my $poly = Math::Primality::BigPolynomial->new( [$a, 1 ]);
+ #print "Start: ", $poly->string, "\n";
+ $poly->powmod($n, $n, $intr);
+ #print "Finish: ", $poly->string, "\n";
+ $compare->setCoef($a, 0);
+ #print "Compare: ", $compare->string, "\n";
- if($res->isEqual($compare)) {
+ if ( ! $poly->isEqual($compare) ) {
debug "Found not prime at $a\n";
return 0;
}
View
208 lib/Math/Primality/BigPolynomial.pm
@@ -17,30 +17,53 @@ sub new {
my $self = {};
my $class = shift;
my $construction_junk = shift;
+
+ $self->{ZERO} = Math::GMPz->new(0);
+
if ($construction_junk) {
my $type = ref $construction_junk;
if ( $type eq 'ARRAY' ) {
- $self->{COEF} = $construction_junk;
- } elsif ( $type eq 'Math::Primality::BigPolynomial') {
- foreach my $coef (@{$construction_junk->{COEF}}) {
- my $temp = Rmpz_init_set($coef);
- push @{$self->{COEF}}, $temp;
- }
+ $self->{COEF} = $construction_junk;
+ } elsif ($type eq 'Math::Primality::BigPolynomial') {
+ copy($self, $construction_junk);
} else {
my $a = [];
for ( my $i = 0 ; $i < $construction_junk ; $i++ ) {
- push @$a, Math::GMPz->new(0);
+ push @$a, $self->{ZERO};
}
$self->{COEF} = $a;
}
}
else {
- $self->{COEF} = [ Math::GMPz->new(0) ];
+ $self->{COEF} = [ $self->{ZERO} ];
}
bless( $self, $class );
return $self;
}
+sub copy {
+ my $self = shift;
+ my $other = shift;
+ die "Copy called with a non-object" unless ref($other) eq 'Math::Primality::BigPolynomial';
+ $self->{COEF} = [ map { Math::GMPz->new($_) } @{$other->{COEF}} ];
+ 1;
+}
+
+sub string {
+ my $self = shift;
+ my @coefs = $self->coef;
+ my $string = '';
+ foreach my $i (reverse 1 .. $#coefs) {
+ my $c = $coefs[$i];
+ next unless $c;
+ $string .= $c if $c != 1;
+ $string .= ($i > 1) ? "x^$i" : "x";
+ $string .= ' + ';
+ }
+ $string .= $coefs[0];
+ return $string;
+}
+
sub coef {
my $self = shift;
if (@_) { @{ $self->{COEF} } = @_ }
@@ -55,23 +78,17 @@ sub degree {
sub getCoef {
my $self = shift;
my $i = shift;
- if ( $i > $self->degree() ) {
- return Math::GMPz->new(0);
- }
- return undef if $i < 0;
+ return if !defined $i || $i < 0;
+ return $self->{ZERO} if $i >= scalar @{$self->{COEF}};
return $self->{COEF}->[$i];
}
sub isEqual {
- my $self = shift;
- my $other_polynomial = shift;
- if ( $self->degree() != $other_polynomial->degree() ) {
- return 0;
- }
- for ( my $i = 0 ; $i < $self->degree() ; $i++ ) {
- if ( $self->getCoef($i) != $other_polynomial->getCoef($i) ) {
- return 0;
- }
+ my $self = shift;
+ my $other = shift;
+ return 0 unless $self->degree == $other->degree;
+ foreach my $i (0 .. $self->degree) {
+ return 0 unless $self->getCoef($i) == $other->getCoef($i);
}
return 1;
}
@@ -80,104 +97,105 @@ sub setCoef {
my $self = shift;
my $new_coef = shift;
my $index = shift;
- if ( $index < 0 ) {
- die "coef is less than 0";
- }
+ die "setCoef: coef not given" unless defined $index;
+ die "setCoef: coef $index is negative" if $index < 0;
- if ( $index > $self->degree() ) {
- for ( my $j = $self->degree() + 1 ; $j < $index ; $j++ ) {
- push @{ $self->{COEF} }, Math::GMPz->new(0);
- }
- $self->{COEF}->[$index] = $new_coef;
- $self->degree($index);
- }
- else {
- $self->{COEF}->[$index] = $new_coef;
+ for ( my $j = $self->degree() + 1 ; $j < $index ; $j++ ) {
+ $self->{COEF}->[$j] = $self->{ZERO};
}
+ $self->{COEF}->[$index] = Math::GMPz->new($new_coef);
}
sub compact {
my $self = shift;
- my $i = 0;
- LOOP: for ( $i = $self->degree(); $i > 0 ; $i-- ) {
- if ( Math::GMPz::Rmpz_cmp_ui( $self->getCoef($i), 0 ) != 0 ) {
- last LOOP;
- }
- pop @{ $self->{COEF} };
- }
- if ( $i != $self->degree() ) {
- $self->degree( $i );
+ while (scalar @{$self->{COEF}} > 1 && $self->{COEF}->[-1] == 0) {
+ pop @{ $self->{COEF} };
}
+ return $self;
}
sub clear {
my $self = shift;
- $self->{COEF} = [ Math::GMPz->new(0) ];
+ $self->{COEF} = [ $self->{ZERO} ];
}
-sub mpz_poly_mod_mult {
- my ( $rop, $copy_x, $copy_y, $mod, $polymod ) = @_;
- my $x = Math::Primality::BigPolynomial->new($copy_x);
- my $y = Math::Primality::BigPolynomial->new($copy_y);
-
- die "mpz_poly_mod_mult: polymod must be defined!" unless $polymod;
-
- $rop->clear();
-
- my $xdeg = ref $x ? $x->degree() : 0;
- my $ydeg = ref $y ? $y->degree() : 0;
+sub mulmod {
+ my ($self, $copy_y, $mod, $polymod ) = @_;
+ die "mulmod: first argument must be a poly!" unless ref($copy_y) eq 'Math::Primality::BigPolynomial';
+ die "mulmod: mod must be defined and > 0!" unless $mod;
+ die "mulmod: polymod must be defined and > 0!" unless $polymod;
+
+ # Bypassing getCoef is a 2-3x speedup.
+ my @x = @{$self->{COEF}};
+ my @y = @{$copy_y->{COEF}};
+ my $xdeg = $#x;
+ my $ydeg = $#y;
my $maxdeg = $xdeg < $ydeg ? $ydeg : $xdeg;
-
- LOOP: for ( my $i = 0 ; $i < $polymod ; $i++ ) {
- my $sum = Math::GMPz->new(0);
- my $temp = Math::GMPz->new(0);
- for ( my $j = 0 ; $j <= $i ; $j++ ) {
- Rmpz_add($temp, $y->getCoef( $i - $j ),
- $y->getCoef( $i + $polymod - $j ) );
- Rmpz_mul( $temp, $x->getCoef($j), $temp );
- Rmpz_add( $sum, $sum, $temp );
+ # Track which coefficients have non-zero components.
+ my @xset = map { Rmpz_sgn($_) > 0 } @x;
+ my @yset = map { Rmpz_sgn($_) > 0 } @y;
+
+ $self->clear();
+
+ # We use these in the loop.
+ my $sum = Math::GMPz->new(0);
+ my $temp = Math::GMPz->new(0);
+
+ for ( my $k = 0 ; $k < $polymod ; $k++ ) {
+ Rmpz_set_ui($sum, 0);
+ for ( my $i = 0 ; $i <= $k ; $i++ ) {
+ # sum += x[i] * (y[k-i] + y[k+polymod-i])
+ # it turns out one, two, or three of the terms above are typically
+ # zero, so only do as much work as needed.
+ next unless $xset[$i];
+ if ($yset[$k-$i]) {
+ if ($yset[$k+$polymod-$i]) {
+ Rmpz_add( $temp, $y[$k-$i], $y[$k+$polymod-$i] );
+ Rmpz_addmul( $sum, $x[$i], $temp );
+ } else {
+ Rmpz_addmul( $sum, $x[$i], $y[$k-$i] );
+ }
+ } else {
+ if ($yset[$k+$polymod-$i]) {
+ Rmpz_addmul( $sum, $x[$i], $y[$k+$polymod-$i] );
+ } else {
+ # Nothing
+ }
+ }
}
- for ( my $j = 0 ; $j < ( $i + $polymod ) ; $j++ ) {
- Rmpz_mul( $temp, $x->getCoef($j),
- $y->getCoef( $i + $polymod - $j ) );
- Rmpz_add( $sum, $sum, $temp );
+ # Only do this loop if there are any non-zero components to multiply.
+ if ( ($k+1) <= $xdeg ) {
+ for ( my $i = $k+1 ; $i <= ( $k + $polymod ) ; $i++ ) {
+ next unless $xset[$i] && $yset[$k+$polymod-$i];
+ Rmpz_addmul( $sum, $x[$i], $y[$k+$polymod-$i] );
+ }
}
Rmpz_mod( $temp, $sum, $mod );
- $rop->setCoef( $temp, $i );
+ $self->setCoef( $temp, $k );
- if ( $i > $maxdeg && Rmpz_cmp_ui( $sum, 0 ) == 0 ) {
- last LOOP;
- }
+ last if $k > $maxdeg && $sum == 0;
}
- $rop->compact();
+ $self->compact();
}
-sub mpz_poly_mod_power {
- my ( $rop, $x, $power, $mult_mod, $poly_mod ) = @_;
-
- die "mpz_poly_mod_power: polymod must be defined!" unless $poly_mod;
-
- $rop->clear();
- $rop->setCoef( Math::GMPz->new(1), 0 );
-
- my $i = Rmpz_sizeinbase( $power, 2 );
-
- LOOP: for ( ; $i >= 0 ; $i-- ) {
- mpz_poly_mod_mult( $rop, $rop, $rop, $mult_mod, $poly_mod );
-
- if ( Rmpz_tstbit( $power, $i ) ) {
- mpz_poly_mod_mult( $rop, $rop, $x, $mult_mod, $poly_mod );
- }
-
- if ( $i == 0 ) {
- last LOOP;
- }
- }
-
- $rop->compact();
+sub powmod {
+ my ($self, $power, $mult_mod, $poly_mod) = @_;
+ die "mpz_poly_mod_power: polymod must be defined!" unless $poly_mod;
+
+ my $x = Math::Primality::BigPolynomial->new($self);
+ $self->clear();
+ $self->setCoef(1, 0);
+
+ my $p = Math::GMPz->new($power);
+ while (Rmpz_sgn($p) > 0) {
+ $self->mulmod($x, $mult_mod, $poly_mod) if Rmpz_odd_p($p);
+ Rmpz_div_2exp($p, $p, 1);
+ $x->mulmod($x, $mult_mod, $poly_mod) if Rmpz_sgn($p) > 0;
+ }
+ $self->compact();
}
1;
Please sign in to comment.
Something went wrong with that request. Please try again.