This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.23
in repository libmath-prime-util-perl.

commit 574e36e8f6f8a42149da5a6a8010f905b73b3968
Author: Dana Jacobsen <d...@acm.org>
Date:   Sat Mar 2 17:16:07 2013 -0800

    Test coverage and small AKS speedup
---
 Changes                |  9 ++++++---
 MANIFEST               |  1 +
 README                 |  2 +-
 XS.xs                  | 16 +++++++++-------
 aks.c                  | 20 ++++++++++++++++----
 lib/Math/Prime/Util.pm |  7 ++++---
 t/19-moebius.t         | 15 +++++++++++++--
 xt/primality-aks.pl    | 28 ++++++++++++++++++++++++++++
 8 files changed, 78 insertions(+), 20 deletions(-)

diff --git a/Changes b/Changes
index 4d3c052..994544b 100644
--- a/Changes
+++ b/Changes
@@ -2,8 +2,6 @@ Revision history for Perl extension Math::Prime::Util.
 
 0.23 xx March 2013
 
-    - exp_mangoldt in XS, and speed up the PP version.
-
     - Replace XS Zeta for x > 5 with series from Cephes.  It is 1 eps more
       accurate for a small fraction of inputs.  More importantly, it is much
       faster in range 5 < x < 10.  This only affects non-integer inputs.
@@ -20,7 +18,12 @@ Revision history for Perl extension Math::Prime::Util.
 
     - Add the first and second Chebyshev functions (theta and psi).
 
-    - Divisor sum with no sub is ~10x faster.
+    - Performance:
+       - Divisor sum with no sub is ~10x faster.
+       - Speed up PP version of exp_mangoldt, create XS version.
+       - Zeta much faster as mentioned above.
+       - faster nth_prime as mentioned above.
+       - AKS about 10% faster.
 
 0.22 26 February 2013
 
diff --git a/MANIFEST b/MANIFEST
index 50ca9e0..8506ef8 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -87,6 +87,7 @@ t/91-release-pod-syntax.t
 t/92-release-pod-coverage.t
 xt/moebius-mertens.pl
 xt/primality-small.pl
+xt/primality-aks.pl
 xt/factor-holf.pl
 xt/make-script-test-data.pl
 .travis.yml
diff --git a/README b/README
index 7ff2d93..9c5a969 100644
--- a/README
+++ b/README
@@ -1,4 +1,4 @@
-Math::Prime::Util version 0.22
+Math::Prime::Util version 0.23
 
 A set of utilities related to prime numbers.  These include multiple sieving
 methods, is_prime, prime_count, nth_prime, approximations and bounds for
diff --git a/XS.xs b/XS.xs
index 490b880..b46185e 100644
--- a/XS.xs
+++ b/XS.xs
@@ -444,11 +444,14 @@ _XS_mertens(IN UV n)
 UV
 _XS_exp_mangoldt(IN UV n)
   CODE:
-    if (n <= 1)           XSRETURN_UV(1);
-    if ((n & (n-1)) == 0) XSRETURN_UV(2);   /* Power of 2 */
-    if ((n & 1) == 0)     XSRETURN_UV(1);   /* Even number */
-    /* if (_XS_is_prime(n))  XSRETURN_UV(n); */
-    {
+    if (n <= 1)
+      RETVAL = 1;
+    else if ((n & (n-1)) == 0)  /* Power of 2 */
+      RETVAL = 2;
+    else if ((n & 1) == 0)      /* Even number */
+      RETVAL = 1;
+    /* else if (_XS_is_prime(n))  RETVAL = n; */
+    else {
       UV factors[MPU_MAX_FACTORS+1];
       UV nfactors, i;
       /* We could try a partial factor, e.g. looking for two small factors */
@@ -458,9 +461,8 @@ _XS_exp_mangoldt(IN UV n)
         if (factors[i] != factors[0])
           XSRETURN_UV(1);
       }
-      XSRETURN_UV(factors[0]);
+      RETVAL = factors[0];
     }
-    RETVAL = 1;
   OUTPUT:
     RETVAL
 
diff --git a/aks.c b/aks.c
index 4e534ce..0e3b123 100644
--- a/aks.c
+++ b/aks.c
@@ -131,16 +131,28 @@ static void poly_mod_mul(UV* px, UV* py, UV* res, UV r, 
UV mod)
 }
 static void poly_mod_sqr(UV* px, UV* res, UV r, UV mod)
 {
-  UV d, s, sum, rindex;
+  UV c, d, s, sum, rindex, maxpx;
   UV degree = r-1;
 
   memset(res, 0, r * sizeof(UV)); /* zero out sums */
+  /* Discover index of last non-zero value in px */
+  for (s = degree; s > 0; s--)
+    if (px[s] != 0)
+      break;
+  maxpx = s;
+  /* 1D convolution */
   for (d = 0; d <= 2*degree; d++) {
+    UV s_beg = (d <= degree) ? 0 : d-degree;
+    UV s_end = ((d/2) <= maxpx) ? d/2 : maxpx;
+    if (s_end < s_beg) continue;
     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];
+    for (s = s_beg; s < s_end; s++) {
+      c = px[s];
+      sum += 2*c * px[d-s];
     }
+    /* Special treatment for last point */
+    c = px[s_end];
+    sum += (s_end*2 == d)  ?  c*c  :  2*c*px[d-s_end];
     rindex = (d < r) ? d : d-r;  /* d % r */
     res[rindex] = (res[rindex] + sum) % mod;
   }
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index e491903..ec4dd9e 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -6,7 +6,7 @@ use Bytes::Random::Secure;
 
 BEGIN {
   $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
-  $Math::Prime::Util::VERSION = '0.22';
+  $Math::Prime::Util::VERSION = '0.23';
 }
 
 # parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
@@ -1045,7 +1045,8 @@ sub consecutive_integer_lcm {
   return 0 if $n < 1;
 
   my $pn = 1;
-  if ($n >= (($_Config{'maxbits'} == 32) ? 22 : 46)) {
+  my $max = ($_Config{'maxbits'} == 32) ? 22 : ($] < 5.008) ? 43 : 46;
+  if ($n >= $max) {
     if (!defined $Math::BigInt::VERSION) {
       eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
       or do { croak "Cannot load Math::BigInt"; };
@@ -1990,7 +1991,7 @@ Math::Prime::Util - Utilities related to prime numbers, 
including fast sieves an
 
 =head1 VERSION
 
-Version 0.22
+Version 0.23
 
 
 =head1 SYNOPSIS
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 433b54e..9b976cd 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -161,12 +161,12 @@ plan tests => 0 + 1
                 + 3*scalar(keys %mertens)
                 + 1*scalar(keys %big_mertens)
                 + 2 # Small Phi
-                + scalar(keys %totients)
+                + 6 + scalar(keys %totients)
                 + scalar(keys %jordan_totients)
                 + 2  # Dedekind psi calculated two ways
                 + 1  # Calculate J5 two different ways
                 + 2 * $use64 # Jordan totient example
-                + scalar(keys %sigmak)
+                + 1 + scalar(keys %sigmak)
                 + scalar(keys %mangoldt)
                 + scalar(keys %chebyshev1)
                 + scalar(keys %chebyshev2);
@@ -204,6 +204,12 @@ while (my($n, $mertens) = each (%big_mertens)) {
 while (my($n, $phi) = each (%totients)) {
   is( euler_phi($n), $phi, "euler_phi($n) == $phi" );
 }
+is_deeply( [euler_phi(0,0)], [0],     "euler_phi(0,0)" );
+is_deeply( [euler_phi(1,0)], [],      "euler_phi with end < start" );
+is_deeply( [euler_phi(0,1)], [0,1],   "euler_phi 0-1" );
+is_deeply( [euler_phi(1,2)], [1,1],   "euler_phi 1-2" );
+is_deeply( [euler_phi(1,3)], [1,1,2], "euler_phi 1-3" );
+is_deeply( [euler_phi(2,3)], [1,2],   "euler_phi 2-3" );
 
 ###### Jordan Totient
 while (my($k, $tref) = each (%jordan_totients)) {
@@ -248,6 +254,11 @@ while (my($k, $sigmaref) = each (%sigmak)) {
   }
   is_deeply( \@slist, $sigmaref, "Sum of divisors to the ${k}th power: 
Sigma_$k" );
 }
+# k=1 standard sum -- much faster
+{
+  my @slist = map { divisor_sum($_) } 1 .. scalar @{$sigmak{1}};
+  is_deeply(\@slist, $sigmak{1}, "divisor_sum(n)");
+}
 
 ###### Exponential of von Mangoldt
 while (my($n, $em) = each (%mangoldt)) {
diff --git a/xt/primality-aks.pl b/xt/primality-aks.pl
new file mode 100755
index 0000000..c0cf7f5
--- /dev/null
+++ b/xt/primality-aks.pl
@@ -0,0 +1,28 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+$| = 1;  # fast pipes
+
+my $limit = shift || 10_000_000;
+
+use Math::Prime::Util qw/is_aks_prime/;
+use Math::Prime::FastSieve;
+my $sieve = Math::Prime::FastSieve::Sieve->new($limit + 10_000);
+
+if (1) {
+  my $n = 2;
+  while ($n <= $limit) {
+    print "$n\n" if $n > 69000; # unless $i++ % 262144;
+    die "$n should be prime" unless is_aks_prime($n);
+    my $next = $sieve->nearest_ge( $n+1 );
+    my $diff = ($next - $n) >> 1;
+    if ($diff > 1) {
+      foreach my $d (1 .. $diff-1) {
+        my $cn = $n + 2*$d;
+        die "$cn should be composite" if is_aks_prime($cn);
+      }
+    }
+    $n = $next;
+  }
+  print "Success to $limit!\n";
+}

-- 
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

Reply via email to