[libmath-prime-util-perl] 12/20: Test coverage and small AKS speedup
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:47:31 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list