[libmath-prime-util-perl] 45/181: Use presieve/test for large segment bases. Much faster and less memory.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:05 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit 0b3812fce8f94a4c36086f724fae40551f31d098
Author: Dana Jacobsen <dana at acm.org>
Date: Wed Dec 25 03:51:51 2013 -0800
Use presieve/test for large segment bases. Much faster and less memory.
---
MANIFEST | 1 +
XS.xs | 7 +------
lib/Math/Prime/Util.pm | 14 +++++++-------
sieve.c | 36 ++++++++++++++++++++----------------
4 files changed, 29 insertions(+), 29 deletions(-)
diff --git a/MANIFEST b/MANIFEST
index 2f77ac4..7f91dd8 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -104,6 +104,7 @@ xt/primality-aks.pl
xt/primality-proofs.pl
xt/small-is-next-prev.pl
xt/factor-holf.pl
+xt/legendre_phi.t
xt/make-script-test-data.pl
xt/measure_zeta_accuracy.pl
xt/pari-totient-moebius.pl
diff --git a/XS.xs b/XS.xs
index 8d17681..bcc417b 100644
--- a/XS.xs
+++ b/XS.xs
@@ -834,12 +834,7 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
* exact limits will change based on the sieve vs. next_prime speed. */
if (beg <= end && !CvISXSUB(cv) && (
#if BITS_PER_WORD == 64
- (beg >= UVCONST(10000000000000000000) && end-beg <1000000000) ||
- (beg >= UVCONST( 1000000000000000000) && end-beg < 95000000) ||
- (beg >= UVCONST( 100000000000000000) && end-beg < 8000000) ||
- (beg >= UVCONST( 10000000000000000) && end-beg < 1700000) ||
- (beg >= UVCONST( 1000000000000000) && end-beg < 400000) ||
- (beg >= UVCONST( 100000000000000) && end-beg < 130000) ||
+ (beg >= UVCONST( 100000000000000) && end-beg < 100000) ||
(beg >= UVCONST( 10000000000000) && end-beg < 40000) ||
(beg >= UVCONST( 1000000000000) && end-beg < 17000) ||
#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index b0ff55d..c4cbf66 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -2857,12 +2857,12 @@ less memory than returning an array directly).
my @primes = @{ primes( 500 ) };
- print "$_\n" for (@{primes( 20, 100 )});
+ print "$_\n" for @{primes(20,100)};
Sieving will be done if required. The algorithm used will depend on the range
-and whether a sieve result already exists. Possibilities include trial
-division (for ranges with only one expected prime), a Sieve of Eratosthenes
-using wheel factorization, or a segmented sieve.
+and whether a sieve result already exists. Possibilities include primality
+testing (for very small ranges), a Sieve of Eratosthenes using wheel
+factorization, or a segmented sieve.
=head2 next_prime
@@ -4866,9 +4866,9 @@ code may be faster for very large values, but isn't available for testing.
Note that the Sieve of Atkin is I<not> faster than the Sieve of Eratosthenes
when both are well implemented. The only Sieve of Atkin that is even
-competitive is Bernstein's super optimized I<primegen>, which runs about
-10% faster than the simple SoE in this module, slower than Pari and yafu's
-SoE implementations, and over 2x slower than primesieve.
+competitive is Bernstein's super optimized I<primegen>, which runs slightly
+slower than the SoE in this module. The SoE's in Pari, yafu, and primesieve
+are all faster.
=item Prime Counts and Nth Prime
diff --git a/sieve.c b/sieve.c
index 75ea093..3f818ca 100644
--- a/sieve.c
+++ b/sieve.c
@@ -9,6 +9,9 @@
#include "cache.h"
#include "util.h"
+/* If the base sieve is larger than this, presieve and test */
+#define BASE_SIEVE_LIMIT 4000000
+
/* 1001 bytes of presieved mod-30 bytes. If the area to be sieved is
* appropriately filled with this data, then 7, 11, and 13 do not have
@@ -284,8 +287,7 @@ unsigned char* sieve_erat30(UV end)
int sieve_segment(unsigned char* mem, UV startd, UV endd)
{
const unsigned char* sieve;
- UV limit;
- UV pcsize;
+ UV limit, slimit;
UV startp = 30*startd;
UV endp = (endd >= (UV_MAX/30)) ? UV_MAX-2 : 30*endd+29;
@@ -298,14 +300,15 @@ int sieve_segment(unsigned char* mem, UV startd, UV endd)
limit = isqrt(endp); /* floor(sqrt(n)), will include p if p*p=endp */
/* Don't use a sieve prime such that p*p > UV_MAX */
if (limit > max_sieve_prime) limit = max_sieve_prime;
- /* printf("segment sieve from %"UVuf" to %"UVuf" (aux sieve to %"UVuf")\n", startp, endp, limit); */
- pcsize = get_prime_cache(limit, &sieve);
- if (pcsize < limit) {
+ slimit = limit;
+ if (slimit > BASE_SIEVE_LIMIT) slimit = BASE_SIEVE_LIMIT;
+ /* printf("segment sieve from %"UVuf" to %"UVuf" (aux sieve to %"UVuf")\n", startp, endp, slimit); */
+ if (get_prime_cache(slimit, &sieve) < slimit) {
release_prime_cache(sieve);
return 0;
}
- START_DO_FOR_EACH_SIEVE_PRIME(sieve, 17, limit)
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, 17, slimit)
{
/* p increments from 17 to at most sqrt(endp). Note on overflow:
* 32-bit: limit= 65535, max p = 65521, p*p = ~0-1965854
@@ -366,8 +369,14 @@ int sieve_segment(unsigned char* mem, UV startd, UV endd)
}
}
END_DO_FOR_EACH_SIEVE_PRIME;
-
release_prime_cache(sieve);
+
+ if (limit > slimit) { /* We've sieved out most composites, but not all. */
+ START_DO_FOR_EACH_SIEVE_PRIME(mem, 0, endp-startp) {
+ if (!_XS_BPSW(startp + p)) /* If the candidate is not prime, */
+ mem[d_] |= mask_; /* mark the sieve location. */
+ } END_DO_FOR_EACH_SIEVE_PRIME;
+ }
return 1;
}
@@ -404,6 +413,7 @@ typedef struct {
void* start_segment_primes(UV low, UV high, unsigned char** segmentmem)
{
segment_context_t* ctx;
+ UV slimit;
MPUassert( high >= low, "start_segment_primes bad arguments");
New(0, ctx, 1, segment_context_t);
@@ -418,17 +428,11 @@ void* start_segment_primes(UV low, UV high, unsigned char** segmentmem)
croak("start_segment_primes: Could not get segment");
*segmentmem = ctx->segment;
-#if 0
- /* If we used this, we would be independent of the primary cache.
- * I think instead, sieve_segment ought to use an aux segmented sieve */
- ctx->base = sieve_erat30( isqrt(ctx->endp)+1 );
- if (ctx->base == 0)
- croak("start_segment_primes: Could not get base");
-#else
ctx->base = 0;
/* Expand primary cache so we won't regen each call */
- get_prime_cache( isqrt(ctx->endp)+1 , 0);
-#endif
+ slimit = isqrt(ctx->endp)+1;
+ if (slimit > BASE_SIEVE_LIMIT) slimit = BASE_SIEVE_LIMIT;
+ get_prime_cache( slimit, 0);
return (void*) ctx;
}
--
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