Permalink
Browse files

WIP clean up AKS stuff

  • Loading branch information...
1 parent f002dec commit 979a82a5e2c8e4462bc70152e81ad6b36eb2deee @leto committed Dec 30, 2011
Showing with 103 additions and 2 deletions.
  1. +2 −0 .gitignore
  2. +98 −0 aks.pl
  3. +2 −2 lib/Math/Primality/BigPolynomial.pm
  4. +1 −0 t/is_aks_prime.t
  5. BIN xt/primes1.zip
View
@@ -23,3 +23,5 @@ xt/*.txt
META.yml
check_dist
nytprof/
+MYMETA.json
+MYMETA.yml
View
@@ -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";
@@ -112,8 +112,8 @@ sub mpz_poly_mod_mult {
$rop->clear();
- my $xdeg = $x->degree();
- my $ydeg = $y->degree();
+ my $xdeg = $x ? $x->degree() : 0;
+ my $ydeg = $y ? $y->degree() : 0;
my $maxdeg = $xdeg < $ydeg ? $ydeg : $xdeg;
LOOP: for ( my $i = 0 ; $i < $polymod ; $i++ ) {
View
@@ -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
Binary file not shown.

0 comments on commit 979a82a

Please sign in to comment.