Browse files

Merge branch 'fix_aks', fixes #1

Conflicts:
	.gitignore
  • Loading branch information...
2 parents a1ea2be + 8a5aca4 commit 1c6545f8f7a464bf5a5066cc4ef0d06532631588 @leto committed Sep 26, 2012
Showing with 191 additions and 53 deletions.
  1. +4 −0 .gitignore
  2. +98 −0 aks.pl
  3. +25 −21 lib/Math/Primality/AKS.pm
  4. +34 −31 lib/Math/Primality/BigPolynomial.pm
  5. BIN nytprof.out
  6. +3 −1 t/00-load.t
  7. +26 −0 t/big_polynomial.t
  8. +1 −0 t/is_aks_prime.t
  9. BIN xt/primes1.zip
View
4 .gitignore
@@ -23,6 +23,10 @@ tags
xt/*.txt
META.yml
check_dist
+<<<<<<< HEAD
MYMETA.json
MYMETA.yml
META.json
+nytprof/
+MYMETA.json
+MYMETA.yml
View
98 aks.pl
@@ -0,0 +1,98 @@
+#!/usr/bin/perl
+
+use strict;
+use warnings;
+
+{
+ package polynomial;
+
+package main;
+
+use Math::GMPz qw/:mpz/;
+use Data::Dumper;
+
+if (scalar(@ARGV) != 1) {
+ print "Usage: aks.pl <number to test>\n";
+ exit
+}
+
+my $n;
+print "Math::GMPz version: $Math::GMPz::VERSION\n";
+eval {
+ $n = Math::GMPz->new($ARGV[0]);
+ 1;
+} or do {
+ print "Error converting ", $ARGV[0], "to a number\n";
+ exit;
+};
+print "Running AKS with number $n\n";
+
+if(Rmpz_perfect_power_p($n)) {
+ print "$n is a perfect power.\n";
+ exit;
+}
+
+my $r = Math::GMPz->new(2);
+my $logn = Rmpz_sizeinbase($n, 2);
+my $limit = Math::GMPz->new($logn * $logn);
+Rmpz_mul_ui($limit, $limit, 4);
+
+# Witness search
+
+OUTERLOOP: while (Rmpz_cmp($r, $n) == -1) {
+ if(Rmpz_divisible_p($n, $r)) {
+ print "$n is divisible by $r\n";
+ exit;
+ }
+
+ if(Rmpz_probab_prime_p($n, 5)) {
+ my $i = Math::GMPz->new(1);
+
+ INNERLOOP: for ( ; Rmpz_cmp($n, $limit) <= 0; Rmpz_add_ui($i, $i, 1)) {
+ my $res = Math::GMPz->new(0);
+ Rmpz_powm($res, $n, $i, $r);
+ if (Rmpz_cmp_ui($res, 1) == 0) {
+ last OUTERLOOP;
+ }
+ }
+
+ }
+ Rmpz_add_ui($r, $r, 1);
+}
+if (Rmpz_cmp($r, $n) == 0) {
+ print "Found $n is prime while checking for r\n";
+ exit;
+}
+
+# Polynomial check
+my $a;
+my $sqrtr = Math::GMPz->new(0);
+
+Rmpz_sqrt($sqrtr, $r);
+my $polylimit = Math::GMPz->new(0);
+Rmpz_add_ui($polylimit, $sqrtr, 1);
+Rmpz_mul_ui($polylimit, $polylimit, $logn);
+Rmpz_mul_ui($polylimit, $polylimit, 2);
+
+my $intr = Rmpz_get_ui($r);
+
+for($a = 1; Rmpz_cmp_ui($polylimit, $a) <= 0; $a++) {
+ print "Checking at $a\n";
+ my $final_size = Math::GMPz->new(0);
+ Rmpz_mod($final_size, $n, $r);
+ my $compare = polynomial->new(Rmpz_get_ui($final_size));
+ $compare->setCoef(1, Rmpz_get_ui($final_size));
+ $compare->setCoef($a, 0);
+ my $res = polynomial->new($intr);
+ my $base = polynomial->new(1);
+ $base->setCoef($a, 0);
+ $base->setCoef(1, 1);
+
+ mpz_poly_mod_power($res, $base, $n, $n, $intr);
+
+ if($res->isEqual($compare)) {
+ print "Found not prime at $a\n";
+ exit;
+ }
+}
+print "Is prime\n";
View
46 lib/Math/Primality/AKS.pm
@@ -3,6 +3,8 @@ use warnings;
use strict;
use Math::GMPz qw/:mpz/;
+use Math::Primality::BigPolynomial;
+
use base 'Exporter';
use Carp qw/croak/;
@@ -35,7 +37,9 @@ our %EXPORT_TAGS = ( all => \@EXPORT_OK );
=head1 FUNCTIONS
-=head2 aks($n)
+=head2 is_aks_prime($n)
+
+Returns 1 if $n is an AKS prime, 0 if it is not.
=cut
@@ -71,13 +75,13 @@ sub is_aks_prime($) {
if(Rmpz_probab_prime_p($n, 5)) {
my $i = Math::GMPz->new(1);
+ my $res = Math::GMPz->new(0);
INNERLOOP: for ( ; Rmpz_cmp($n, $limit) <= 0; Rmpz_add_ui($i, $i, 1)) {
- my $res = Math::GMPz->new(0);
- Rmpz_powm($res, $n, $i, $r);
- if (Rmpz_cmp_ui($res, 1) == 0) {
- last OUTERLOOP;
- }
+ Rmpz_powm($res, $n, $i, $r);
+ if (Rmpz_cmp_ui($res, 1) == 0) {
+ last OUTERLOOP;
+ }
}
}
@@ -90,7 +94,7 @@ sub is_aks_prime($) {
# Polynomial check
my $a;
- my $sqrtr = Math::GMPz->new(0);
+ my $sqrtr = Math::GMPz->new($r);
Rmpz_sqrt($sqrtr, $r);
my $polylimit = Math::GMPz->new(0);
@@ -100,19 +104,20 @@ sub is_aks_prime($) {
my $intr = Rmpz_get_ui($r);
- for($a = 1; Rmpz_cmp_ui($polylimit, $a) <= 0; $a++) {
+ for($a = 1; Rmpz_cmp_ui($polylimit, $a) >= 0; $a++) {
debug "Checking at $a\n";
my $final_size = Math::GMPz->new(0);
Rmpz_mod($final_size, $n, $r);
- my $compare = polynomial->new(Rmpz_get_ui($final_size));
- $compare->setCoef(1, Rmpz_get_ui($final_size));
- $compare->setCoef($a, 0);
- my $res = polynomial->new($intr);
- my $base = polynomial->new(1);
- $base->setCoef($a, 0);
- $base->setCoef(1, 1);
+ 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);
- mpz_poly_mod_power($res, $base, $n, $n, $intr);
if($res->isEqual($compare)) {
debug "Found not prime at $a\n";
@@ -129,11 +134,10 @@ Jonathan "Duke" Leto C<< <jonathan@leto.net> >>
=head1 BUGS
-Please report any bugs or feature requests to C<bug-math-primality-aks at rt.cpan.org>,
-or through the web interface at
-L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Math::Primality::AKS>. I will be
-notified, and then you'll automatically be notified of progress on your bug as I
-make changes.
+Please report any bugs or feature requests to
+C<bug-math-primality-aks at rt.cpan.org>, or through the web interface at
+L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Math::Primality::AKS>.
+I will be notified, and then you'll automatically be notified of progress on your bug as I make changes.
=head1 THANKS
View
65 lib/Math/Primality/BigPolynomial.pm
@@ -2,7 +2,7 @@ package Math::Primality::BigPolynomial;
use strict;
use warnings;
-use Math::GMPz;
+use Math::GMPz qw/:mpz/;
# ABSTRACT: Big Polynomials
@@ -18,22 +18,24 @@ sub new {
my $class = shift;
my $construction_junk = shift;
if ($construction_junk) {
- if ( ref($construction_junk) eq 'ARRAY' ) {
+ my $type = ref $construction_junk;
+ if ( $type eq 'ARRAY' ) {
$self->{COEF} = $construction_junk;
- $self->{DEGREE} = scalar(@$construction_junk);
- }
- else {
- $self->{DEGREE} = $construction_junk;
- my @a = [];
+ } elsif ( $type eq 'Math::Primality::BigPolynomial') {
+ foreach my $coef (@{$construction_junk->{COEF}}) {
+ my $temp = Rmpz_init_set($coef);
+ push @{$self->{COEF}}, $temp;
+ }
+ } else {
+ my $a = [];
for ( my $i = 0 ; $i < $construction_junk ; $i++ ) {
- push @a, Math::GMPz->new(0);
+ push @$a, Math::GMPz->new(0);
}
- $self->{COEF} = @a;
+ $self->{COEF} = $a;
}
}
else {
$self->{COEF} = [ Math::GMPz->new(0) ];
- $self->{DEGREE} = 1;
}
bless( $self, $class );
return $self;
@@ -47,17 +49,17 @@ sub coef {
sub degree {
my $self = shift;
- if (@_) { $self->{DEGREE} = shift }
- return $self->{DEGREE};
+ return (scalar @{$self->{COEF}} - 1);
}
sub getCoef {
my $self = shift;
my $i = shift;
if ( $i > $self->degree() ) {
- return 0;
+ return Math::GMPz->new(0);
}
- return ${ $self->{COEF} }[$i];
+ return undef if $i < 0;
+ return $self->{COEF}->[$i];
}
sub isEqual {
@@ -83,55 +85,55 @@ sub setCoef {
}
if ( $index > $self->degree() ) {
- for ( my $j = $self->degree() + 1 ; $j <= $index ; $j++ ) {
+ for ( my $j = $self->degree() + 1 ; $j < $index ; $j++ ) {
push @{ $self->{COEF} }, Math::GMPz->new(0);
}
- push @{ $self->{COEF} }, $new_coef;
+ $self->{COEF}->[$index] = $new_coef;
$self->degree($index);
}
else {
- ${ $self->{COEF} }[$index] = $new_coef;
+ $self->{COEF}->[$index] = $new_coef;
}
}
sub compact {
my $self = shift;
my $i = 0;
- LOOP: for ( $i = $self->degree() - 1 ; $i > 0 ; $i-- ) {
+ 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 + 1 );
+ $self->degree( $i );
}
}
sub clear {
my $self = shift;
- $self->{COEF} = undef;
- $self->{DEGREE} = undef;
$self->{COEF} = [ Math::GMPz->new(0) ];
- $self->{DEGREE} = 1;
}
sub mpz_poly_mod_mult {
- my $self = shift;
- my ( $rop, $x, $y, $mod, $polymod ) = @_;
+ 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 = $x->degree();
- my $ydeg = $y->degree();
+ my $xdeg = ref $x ? $x->degree() : 0;
+ my $ydeg = ref $y ? $y->degree() : 0;
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_add($temp, $y->getCoef( $i - $j ),
+ $y->getCoef( $i + $polymod - $j ) );
Rmpz_mul( $temp, $x->getCoef($j), $temp );
Rmpz_add( $sum, $sum, $temp );
}
@@ -154,19 +156,20 @@ sub mpz_poly_mod_mult {
}
sub mpz_poly_mod_power {
- my $self = shift;
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-- ) {
- Rmpz_poly_mod_mult( $rop, $rop, $rop, $mult_mod, $poly_mod );
+ mpz_poly_mod_mult( $rop, $rop, $rop, $mult_mod, $poly_mod );
if ( Rmpz_tstbit( $power, $i ) ) {
- mpz_poly_mod_mul( $rop, $rop, $x, $mult_mod, $poly_mod );
+ mpz_poly_mod_mult( $rop, $rop, $x, $mult_mod, $poly_mod );
}
if ( $i == 0 ) {
View
BIN nytprof.out
Binary file not shown.
View
4 t/00-load.t
@@ -1,10 +1,12 @@
#!/usr/bin/env perl
use strict;
-use Test::More tests => 1;
+use Test::More tests => 3;
BEGIN {
use_ok( 'Math::Primality' );
+ use_ok( 'Math::Primality::AKS' );
+ use_ok( 'Math::Primality::BigPolynomial' );
}
diag( "Testing Math::Primality $Math::Primality::VERSION, Perl $], $^X" );
View
26 t/big_polynomial.t
@@ -0,0 +1,26 @@
+#!/usr/bin/evn perl
+
+use strict;
+use warnings;
+use Test::More;
+#use Carp::Always;
+
+use Math::GMPz qw/:mpz/;
+
+BEGIN {
+ use_ok ('Math::Primality::BigPolynomial' );
+}
+
+
+my $b = Math::Primality::BigPolynomial->new([1,3,7]);
+isa_ok($b, 'Math::Primality::BigPolynomial');
+
+is($b->getCoef(0), 1, 'coef(0) is 1');
+is($b->getCoef(1), 3, 'coef(1) is 3');
+is($b->getCoef(2), 7, 'coef(2) is 7');
+# this is a bit wonky
+cmp_ok($b->getCoef(3),'==', Math::GMPz->new(0), 'coef(3) is 0');
+is($b->getCoef(-1), undef, 'coef(-1) is undef');
+is($b->getCoef(-10), undef, 'coef(-10) is undef');
+
+done_testing;
View
1 t/is_aks_prime.t
@@ -3,6 +3,7 @@
use strict;
use warnings;
use Test::More tests => 232;
+#use Carp::Always;
use Math::GMPz qw/:mpz/;
use POSIX qw(ceil floor); # for testing _Rmpz_logbase2* functions
View
BIN xt/primes1.zip
Binary file not shown.

0 comments on commit 1c6545f

Please sign in to comment.