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 <[email protected]> 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 [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
