This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.13
in repository libmath-prime-util-perl.
commit e440b5da107e790f0c9876b26257368378d687aa
Author: Dana Jacobsen <d...@acm.org>
Date: Sat Nov 17 21:53:56 2012 -0800
Add AKS primality test
---
.gitignore | 2 +
Changes | 6 +-
MANIFEST | 2 +
Makefile.PL | 1 +
TODO | 6 +-
XS.xs | 4 +
aks.c | 210 ++++++++++++++++++++++++++++++++++++++++++++++
aks.h | 8 ++
lib/Math/Prime/Util.pm | 43 ++++++++--
lib/Math/Prime/Util/PP.pm | 156 ++++++++++++++++++++++++++++++++++
sieve.c | 2 +-
t/10-isprime.t | 11 ++-
12 files changed, 437 insertions(+), 14 deletions(-)
diff --git a/.gitignore b/.gitignore
index 4ed058c..04fde6f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,10 +3,12 @@ Makefile
Makefile.old
MYMETA.*
Util.bs
+pm_to_blib
XS.[co]
cache.o
factor.o
sieve.o
util.o
lehmer.o
+aks.o
blib*
diff --git a/Changes b/Changes
index 64480bb..3a29363 100644
--- a/Changes
+++ b/Changes
@@ -1,6 +1,6 @@
Revision history for Perl extension Math::Prime::Util.
-0.12 12 November 2012
+0.12 17 November 2012
- Add bin/primes.pl and bin/factor.pl
@@ -10,6 +10,7 @@ Revision history for Perl extension Math::Prime::Util.
prime_set_config set config options
RiemannZeta export and make accurate for small reals
is_provable_prime prove primes after BPSW
+ is_aks_prime prove prime via AKS
- Add 'assume_rh' configuration option (default: false) which can be set
to allow functions to assume the Riemann Hypothesis.
@@ -30,7 +31,8 @@ Revision history for Perl extension Math::Prime::Util.
- bug fix for prime_count on edge of cache.
- - Add Lehmer prime counting algorithm. This is much faster than sieving.
+ - prime_count will use Lehmer prime counting algorithm for largish
+ sizes (above 4 million). This is MUCH faster than sieving.
0.11 23 July 2012
- Turn off threading tests on Cygwin, as threads on some Cygwin platforms
diff --git a/MANIFEST b/MANIFEST
index 608078e..bcc1329 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -10,6 +10,8 @@ README
TODO
XS.xs
ptypes.h
+aks.h
+aks.c
cache.h
cache.c
factor.h
diff --git a/Makefile.PL b/Makefile.PL
index ca13682..ad34e57 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -9,6 +9,7 @@ WriteMakefile1(
OBJECT => 'cache.o ' .
'factor.o ' .
+ 'aks.o ' .
'lehmer.o ' .
'sieve.o ' .
'util.o ' .
diff --git a/TODO b/TODO
index d20aa19..7da582a 100644
--- a/TODO
+++ b/TODO
@@ -23,8 +23,6 @@
- tests for primorial
-- document prime_set_config
-
- Test all routines for numbers on word-size boundary, or ranges that cross.
- Test all functions return either native or bigints. Functions that return
@@ -32,4 +30,6 @@
- Make proper pminus1 in PP code, like factor.c.
-- Add Lehmer prime count function for large values.
+- Add Lehmer in PP (I have a version, just needs some polishing).
+
+- Move AKS tests into their own test, and add some provable prime tests.
diff --git a/XS.xs b/XS.xs
index ba7d0b3..d61f012 100644
--- a/XS.xs
+++ b/XS.xs
@@ -9,6 +9,7 @@
#include "util.h"
#include "factor.h"
#include "lehmer.h"
+#include "aks.h"
MODULE = Math::Prime::Util PACKAGE = Math::Prime::Util
@@ -49,6 +50,9 @@ _XS_nth_prime(IN UV n)
int
_XS_is_prime(IN UV n)
+int
+_XS_is_aks_prime(IN UV n)
+
UV
_XS_next_prime(IN UV n)
diff --git a/aks.c b/aks.c
new file mode 100644
index 0000000..2b43943
--- /dev/null
+++ b/aks.c
@@ -0,0 +1,210 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+/*
+ * The AKS v6 algorithm, for native integers, where n < 2^(wordbits/2)-1.
+ * Hence on 64-bit machines this works for n < 4294967295, because we do
+ * r = (r + a * b) % n
+ * where r, a, and b are mod n. This could be extended to a full word by
+ * using a mulmod function (like factor.c has), but it's easier to go to
+ * GMP at that point, which also lets one do r or 2r modulos instead of r*r.
+ *
+ * Copyright 2012, Dana Jacobsen.
+ */
+
+#define SQRTN_SHORTCUT 1
+#define VERBOSE 0
+
+#include "ptypes.h"
+#include "util.h"
+#include "sieve.h"
+#include "factor.h"
+#include "cache.h"
+
+#define addmod(n,a,m) ((((m)-(n)) > (a)) ? ((n)+(a)) : ((n)+(a)-(m)))
+
+static UV log2floor(UV n) {
+ UV log2n = 0;
+ while (n >>= 1)
+ log2n++;
+ return log2n;
+}
+
+/* See Bach and Sorenson (1993) for much better algorithm */
+static int is_perfect_power(UV x) {
+ UV b, last;
+ if ((x & (x-1)) == 0) return 1; /* powers of 2 */
+ b = sqrt(x); if (b*b == x) return 1; /* perfect square */
+ b = cbrt(x); if (b*b*b == x) return 1; /* perfect cube */
+ last = log2floor(x) + 1;
+ for (b = 5; b < last; b = _XS_next_prime(b)) {
+ UV root = pow(x, 1.0 / (double)b);
+ if (pow(root, b) == x) return 1;
+ }
+ return 0;
+}
+
+static UV order(UV r, UV n, UV limit) {
+ UV j;
+ UV t = 1;
+ for (j = 1; j <= limit; j++) {
+ t = (t * n) % r;
+ if (t == 1)
+ break;
+ }
+ return j;
+}
+
+static void poly_print(UV* poly, UV r)
+{
+ int i;
+ for (i = r-1; i >= 1; i--) {
+ if (poly[i] != 0)
+ printf("%lux^%d + ", poly[i], i);
+ }
+ if (poly[0] != 0) printf("%lu", poly[0]);
+ printf("\n");
+}
+
+static void poly_mod_mul(UV* px, UV* py, UV* res, UV r, UV mod)
+{
+ int i, j;
+ UV pxi, pyj, rindex;
+
+ memset(res, 0, r * sizeof(UV));
+ for (i = 0; i < r; i++) {
+ pxi = px[i];
+ if (pxi == 0) continue;
+ for (j = 0; j < r; j++) {
+ pyj = py[j];
+ if (pyj == 0) continue;
+ rindex = (i+j) < r ? i+j : i+j-r; /* (i+j) % r */
+ res[rindex] = (res[rindex] + (pxi*pyj) ) % mod;
+ }
+ }
+ memcpy(px, res, r * sizeof(UV)); /* put result in px */
+}
+static void poly_mod_sqr(UV* px, UV* res, UV r, UV mod)
+{
+ int d, s;
+ UV sum, rindex;
+ UV degree = r-1;
+
+ /* we sum a max of r*mod*mod between modulos */
+ if (mod > sqrt(UV_MAX/r))
+ return poly_mod_mul(px, px, res, r, mod);
+
+ memset(res, 0, r * sizeof(UV)); /* zero out sums */
+ for (d = 0; d <= 2*degree; d++) {
+ sum = 0;
+ for (s = (d <= degree) ? 0 : d-degree; s <= (d/2); s++) {
+ UV c = px[s];
+ sum += (s*2 == d) ? c*c : 2*c * px[d-s];
+ }
+ rindex = (d < r) ? d : d-r; /* d % r */
+ res[rindex] = (res[rindex] + sum) % mod;
+ }
+ memcpy(px, res, r * sizeof(UV)); /* put result in px */
+}
+
+static UV* poly_mod_pow(UV* pn, UV power, UV r, UV mod)
+{
+ UV* res;
+ UV* temp;
+
+ Newz(0, res, r, UV);
+ New(0, temp, r, UV);
+ if ( (res == 0) || (temp == 0) )
+ croak("Couldn't allocate space for polynomial of degree %lu\n", r);
+
+ res[0] = 1;
+
+ while (power) {
+ if (power & 1) poly_mod_mul(res, pn, temp, r, mod);
+ power >>= 1;
+ if (power) poly_mod_sqr(pn, temp, r, mod);
+ }
+ Safefree(temp);
+ return res;
+}
+
+static int test_anr(UV a, UV n, UV r)
+{
+ UV* pn;
+ UV* res;
+ int i;
+ int retval = 1;
+
+ Newz(0, pn, r, UV);
+ if (pn == 0)
+ croak("Couldn't allocate space for polynomial of degree %lu\n", r);
+ a %= r;
+ pn[0] = a;
+ pn[1] = 1;
+ res = poly_mod_pow(pn, n, r, n);
+ res[n % r] = addmod(res[n % r], n - 1, n);
+ res[0] = addmod(res[0], n - a, n);
+
+ for (i = 0; i < r; i++)
+ if (res[i] != 0)
+ retval = 0;
+ Safefree(res);
+ Safefree(pn);
+ return retval;
+}
+
+int _XS_is_aks_prime(UV n)
+{
+ UV sqrtn, limit, r, rlimit, a;
+ double log2n;
+
+ /* Check for overflow */
+#if BITS_PER_WORD == 32
+ if (n >= UVCONST(65535))
+#else
+ if (n >= UVCONST(4294967295))
+#endif
+ croak("aks(%"UVuf") overflow", n);
+
+ if (n < 2)
+ return 0;
+ if (n == 2)
+ return 1;
+
+ if (is_perfect_power(n))
+ return 0;
+
+ sqrtn = sqrt(n);
+ log2n = log(n) / log(2); /* C99 has a log2() function */
+ limit = (UV) floor(log2n * log2n);
+
+ if (VERBOSE) { printf("limit is %lu\n", limit); }
+
+ for (r = 2; r < n; r++) {
+ if ((n % r) == 0)
+ return 0;
+#if SQRTN_SHORTCUT
+ if (r > sqrtn)
+ return 1;
+#endif
+ if (order(r, n, limit) > limit)
+ break;
+ }
+
+ if (r >= n)
+ return 1;
+
+ rlimit = (UV) floor(sqrt(r-1) * log2n);
+
+ if (VERBOSE) { printf("r = %lu rlimit = %lu\n", r, rlimit); }
+
+ for (a = 1; a <= rlimit; a++) {
+ if (! test_anr(a, n, r) )
+ return 0;
+ if (VERBOSE) { printf("."); fflush(stdout); }
+ }
+ if (VERBOSE) { printf("\n"); }
+ return 1;
+}
diff --git a/aks.h b/aks.h
new file mode 100644
index 0000000..e8c7452
--- /dev/null
+++ b/aks.h
@@ -0,0 +1,8 @@
+#ifndef MPU_AKS_H
+#define MPU_AKS_H
+
+#include "ptypes.h"
+
+extern int _XS_is_aks_prime(UV n);
+
+#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index f94cbef..eafa8cc 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -16,6 +16,7 @@ our @EXPORT_OK = qw(
prime_precalc prime_memfree
is_prime is_prob_prime is_provable_prime
is_strong_pseudoprime is_strong_lucas_pseudoprime
+ is_aks_prime
miller_rabin
primes
next_prime prev_prime
@@ -819,6 +820,19 @@ sub is_prime {
return is_prob_prime($n);
}
+sub is_aks_prime {
+ my($n) = @_;
+ return 0 if $n <= 0;
+ _validate_positive_integer($n);
+
+ my $max_xs = ($_Config{'maxparam'} > 4294967295) ? 4294967294 : 65534;
+ return _XS_is_aks_prime($n) if $n <= $_XS_MAXVAL && $n <= $max_xs;
+ return Math::Prime::Util::GMP::is_aks_prime($n) if $_HAVE_GMP
+ && defined &Math::Prime::Util::GMP::is_aks_prime;
+ return Math::Prime::Util::PP::is_aks_prime($n);
+}
+
+
sub next_prime {
my($n) = @_;
_validate_positive_integer($n);
@@ -1017,7 +1031,7 @@ sub is_provable_prime {
# prove it. We should do ECPP or APR-CL now, or failing that, do the
# Brillhart-Lehmer-Selfridge test, or Pocklington-Lehmer. Until those
# are written here, we'll do a Lucas test, which is super simple but may
- # be very slow. We could do AKS, but that's even slower.
+ # be very slow. We have AKS code, but it's insanely slow.
# See http://primes.utm.edu/prove/merged.html or other sources.
# It shouldn't be possible to get here without BigInt already loaded.
@@ -1662,10 +1676,10 @@ base under C<10^14>.
Lehmer's method has complexity approximately C<O(b^0.7) + O(a^0.7)>. It
does use more memory however. A calculation of C<Pi(10^14)> completes in
-under 1 minute, C<Pi(10^15)> in approximately 5 minutes, and C<Pi(10^16)>
-in about 30 minutes, however using nearly 800MB of peak memory for the last.
-In contrast, even primesieve using multiple cores would take over a week
-on this same computer to determine C<Pi(10^16)>.
+under 1 minute, C<Pi(10^15)> in undef 5 minutes, and C<Pi(10^16)> in under
+30 minutes, however using nearly 1400MB of peak memory for the last.
+In contrast, even primesieve using 12 cores would take over a week on this
+same computer to determine C<Pi(10^16)>.
Also see the function L</"prime_count_approx"> which gives a very good
approximation to the prime count, and L</"prime_count_lower"> and
@@ -1844,6 +1858,21 @@ usually take less time, but can still fail. Hence you
should always test that
the result is C<2> to ensure the prime is proven.
+=head2 is_aks_prime
+
+ say "$n is definitely prime" if is_aks_prime($n);
+
+Takes a positive number as input, and returns 1 if the input passes the
+Agrawal-Kayal-Saxena (AKS) primality test. This is a deterministic
+unconditional primality test which runs in polynomial time for general input.
+
+This function is only included for completeness and as an example. While the
+implementation is fast compared to the only other Perl implementation available
+(in L<Math::Primality>), it is slow compared to others. However, even
+optimized AKS implementations are far slower than ECPP or other modern
+primality tests.
+
+
=head2 moebius
say "$n is square free" if moebius($n) != 0;
@@ -2293,7 +2322,9 @@ Perl threads.
=head1 PERFORMANCE
Counting the primes to C<10^10> (10 billion), with time in seconds.
-Pi(10^10) = 455,052,511.
+Pi(10^10) = 455,052,511. Note that the Lehmer prime counting method in
+Math::Prime::Util version 0.12 takes 0.064 seconds to do this -- the numbers
+below are all for doing it the long way (sieving).
External C programs in C / C++:
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 181f316..0a67888 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -533,6 +533,36 @@ sub _gcd_ui {
$x;
}
+sub _is_perfect_power {
+ my $n = shift;
+ my $log2n = _log2($n);
+ $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
+ for my $e (@{primes($log2n)}) {
+ return 1 if $n->copy()->broot($e)->bpow($e) == $n;
+ }
+ 0;
+}
+
+sub _order {
+ my($r, $n, $lim) = @_;
+ $lim = $r unless defined $lim;
+
+ return 1 if ($n % $r) == 1;
+ for (my $j = 2; $j <= $lim; $j++) {
+ return $j if _powmod($n, $j, $r) == 1;
+ }
+ return $lim+1;
+}
+
+# same result as: int($n->blog(2)->floor->bstr)
+sub _log2 {
+ my $n = shift;
+ my $log2n = 0;
+ $log2n++ while ($n >>= 1);
+ $log2n;
+}
+
+
sub miller_rabin {
my($n, @bases) = @_;
@@ -783,6 +813,126 @@ sub is_strong_lucas_pseudoprime {
}
+my $_poly_bignum;
+sub _poly_new {
+ my @poly;
+ if ($_poly_bignum) {
+ foreach my $c (@_) {
+ push @poly, (ref $c eq 'Math::BigInt') ? $c->copy :
Math::BigInt->new("$c");
+ }
+ } else {
+ push @poly, $_ for (@_);
+ push @poly, 0 unless scalar @poly;
+ }
+ return \@poly;
+}
+
+sub _poly_print {
+ my($poly) = @_;
+ warn "poly has null top degree" if $#$poly > 0 && !$poly->[-1];
+ foreach my $d (reverse 1 .. $#$poly) {
+ my $coef = $poly->[$d];
+ print "", ($coef != 1) ? $coef : "", ($d > 1) ? "x^$d" : "x", " + "
+ if $coef;
+ }
+ my $p0 = $poly->[0] || 0;
+ print "$p0\n";
+}
+
+sub _poly_mod_mul {
+ my($px, $py, $r, $n) = @_;
+
+ my $px_degree = $#$px;
+ my $py_degree = $#$py;
+ my @res;
+
+ # convolve(px, py) mod (X^r-1,n)
+ my @indices_y = grep { $py->[$_] } (0 .. $py_degree);
+ for (my $ix = 0; $ix <= $px_degree; $ix++) {
+ my $px_at_ix = $px->[$ix];
+ next unless $px_at_ix;
+ foreach my $iy (@indices_y) {
+ my $py_at_iy = $py->[$iy];
+ my $rindex = ($ix + $iy) % $r; # reduce mod X^r-1
+ if (!defined $res[$rindex]) {
+ $res[$rindex] = $_poly_bignum ? Math::BigInt->bzero : 0
+ }
+ $res[$rindex] = ($res[$rindex] + ($py_at_iy * $px_at_ix)) % $n;
+ }
+ }
+ # In case we had upper terms go to zero after modulo, reduce the degree.
+ pop @res while !$res[-1];
+ return \@res;
+}
+
+sub _poly_mod_pow {
+ my($pn, $power, $r, $mod) = @_;
+ my $res = _poly_new(1);
+ my $p = $power;
+
+ while ($p) {
+ $res = _poly_mod_mul($res, $pn, $r, $mod) if ($p & 1);
+ $p >>= 1;
+ $pn = _poly_mod_mul($pn, $pn, $r, $mod) if $p;
+ }
+ return $res;
+}
+
+sub _test_anr {
+ my($a, $n, $r) = @_;
+ my $pp = _poly_mod_pow(_poly_new($a, 1), $n, $r, $n);
+ $pp->[$n % $r] = (($pp->[$n % $r] || 0) - 1) % $n; # subtract X^(n%r)
+ $pp->[ 0] = (($pp->[ 0] || 0) - $a) % $n; # subtract a
+ return 0 if scalar grep { $_ } @$pp;
+ 1;
+}
+
+sub is_aks_prime {
+ my $n = shift;
+ $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
+
+ return 0 if $n < 2;
+ return 0 if _is_perfect_power($n);
+
+ # limit = floor( log2(n) * log2(n) ). o_r(n) must be larger than this
+ my $sqrtn = int(Math::BigFloat->new($n)->bsqrt->bfloor->bstr);
+ my $log2n = Math::BigFloat->new($n)->blog(2);
+ my $log2_squared_n = $log2n * $log2n;
+ my $limit = int($log2_squared_n->bfloor->bstr);
+
+ my $r = next_prime($limit);
+ foreach my $f (@{primes(0,$r-1)}) {
+ return 1 if $f == $n;
+ return 0 if !($n % $f);
+ }
+
+ while ($r < $n) {
+ return 0 if !($n % $r);
+ #return 1 if $r >= $sqrtn;
+ last if _order($r, $n, $limit) > $limit;
+ $r = next_prime($r);
+ }
+
+ return 1 if $r >= $n;
+
+ # Since r is a prime, phi(r) = r-1
+ my $rlimit = int( Math::BigFloat->new($r)->bsub(1)
+ ->bsqrt->bmul($log2n)->bfloor->bstr);
+
+ $_poly_bignum = 1;
+ if ( $n < ( (~0 == 4294967295) ? 65535 : 4294967295 ) ) {
+ $_poly_bignum = 0;
+ $n = int($n->bstr) if ref($n) eq 'Math::BigInt';
+ }
+
+ for (my $a = 1; $a <= $rlimit; $a++) {
+ return 0 unless _test_anr($a, $n, $r);
+ }
+
+ return 1;
+}
+
+
sub _basic_factor {
# MODIFIES INPUT SCALAR
return ($_[0]) if $_[0] < 4;
@@ -1539,6 +1689,12 @@ sources call this a strong Lucas-Selfridge pseudoprime).
This is one half
of the BPSW primality test (the Miller-Rabin strong pseudoprime test with
base 2 being the other half).
+=head2 is_aks_prime
+
+Takes a positive number as input, and returns 1 if the input can be proven
+prime using the AKS primality test. This code is included for completeness
+and as an example, but is impractically slow.
+
=head1 UTILITY FUNCTIONS
diff --git a/sieve.c b/sieve.c
index 138ca24..eb38b82 100644
--- a/sieve.c
+++ b/sieve.c
@@ -186,7 +186,7 @@ int sieve_segment(unsigned char* mem, UV startd, UV endd)
MPUassert( (mem != 0) && (endd >= startd) && (endp >= startp),
"sieve_segment bad arguments");
-
+
/* Fill buffer with marked 7, 11, and 13 */
sieve_prefill(mem, startd, endd);
diff --git a/t/10-isprime.t b/t/10-isprime.t
index a1df315..cdddc9c 100644
--- a/t/10-isprime.t
+++ b/t/10-isprime.t
@@ -3,13 +3,13 @@ use strict;
use warnings;
use Test::More;
-use Math::Prime::Util qw/is_prime/;
+use Math::Prime::Util qw/is_prime is_aks_prime/;
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
my $broken64 = (18446744073709550592 == ~0);
-plan tests => 6 + 19 + 3573 + (5 + 29 + 22 + 23 + 16) + 15 + 27 + 1
+plan tests => 6 + 19 + 3573 + (5 + 29 + 22 + 23 + 16) + 15 + 27 + 1 + 3
+ ($use64 ? 5+1 : 0)
+ ($extra ? 6 : 0)
+ (($extra && $use64) ? 19 : 0);
@@ -118,3 +118,10 @@ SKIP: {
eval { is_prime( $use64 ? "18446744073709551629" : "4294967306" ); };
like($@, qr/range/i, "is_prime on ~0 + delta without bigint should croak");
}
+
+# Simple AKS
+is( is_aks_prime(877), 1, "is_aks_prime(877) is true" );
+# The first number that makes it past the sqrt test to actually run.
+is( is_aks_prime(69197), 1, "is_aks_prime(69197) is true" );
+# A composite
+is( is_aks_prime(1262952907), 0, "is_aks_prime(1262952907) is false" );
--
Alioth's /usr/local/bin/git-commit-notice on
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
