[libmath-prime-util-perl] 02/15: segmented sieve infrastructure, and have forprimes use it.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:48:46 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.29
in repository libmath-prime-util-perl.
commit ede3138f8237707b307356ab9cc3a28685446c95
Author: Dana Jacobsen <dana at acm.org>
Date: Mon May 27 00:02:20 2013 -0700
segmented sieve infrastructure, and have forprimes use it.
---
XS.xs | 38 +++++++++++++++-------
sieve.c | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sieve.h | 3 ++
t/32-iterators.t | 28 +++++++++++++++-
util.c | 2 ++
5 files changed, 155 insertions(+), 13 deletions(-)
diff --git a/XS.xs b/XS.xs
index 3422caa..9bc19c5 100644
--- a/XS.xs
+++ b/XS.xs
@@ -247,6 +247,7 @@ segment_primes(IN UV low, IN UV high);
PREINIT:
AV* av = newAV();
CODE:
+ /* Could rewrite using {start/next/end}_segment_primes functions */
if ((low <= 2) && (high >= 2)) { av_push(av, newSVuv( 2 )); }
if ((low <= 3) && (high >= 3)) { av_push(av, newSVuv( 3 )); }
if ((low <= 5) && (high >= 5)) { av_push(av, newSVuv( 5 )); }
@@ -650,6 +651,9 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
HV *stash;
CV *cv;
SV* svarg;
+ void* ctx;
+ unsigned char* segment;
+ UV seg_base, seg_low, seg_high;
if (!_validate_int(svbeg, 0) || (items >= 3 && !_validate_int(svend,0))) {
dSP;
@@ -675,29 +679,39 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
croak("Not a subroutine reference");
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
+ ctx = start_segment_primes(beg, end, &segment);
if (!CvISXSUB(cv)) {
dMULTICALL;
I32 gimme = G_VOID;
PUSH_MULTICALL(cv);
- START_DO_FOR_EACH_PRIME(beg, end) {
- sv_setuv(svarg, p);
- GvSV(PL_defgv) = svarg;
- MULTICALL;
- } END_DO_FOR_EACH_PRIME
+ while (beg < 7) {
+ beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
+ { sv_setuv(svarg, beg); GvSV(PL_defgv) = svarg; MULTICALL; }
+ beg += 1 + (beg > 2);
+ }
+ while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
+ START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base )
+ { sv_setuv(svarg, seg_base + p); GvSV(PL_defgv) = svarg; MULTICALL; }
+ END_DO_FOR_EACH_SIEVE_PRIME
+ }
#ifdef PERL_HAS_BAD_MULTICALL_REFCOUNT
if (CvDEPTH(multicall_cv) > 1)
SvREFCNT_inc(multicall_cv);
#endif
POP_MULTICALL;
} else {
- START_DO_FOR_EACH_PRIME(beg, end) {
- dSP;
- sv_setuv(svarg, p);
- GvSV(PL_defgv) = svarg;
- PUSHMARK(SP);
- call_sv((SV*)cv, G_VOID|G_DISCARD);
- } END_DO_FOR_EACH_PRIME
+ while (beg < 7) {
+ beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
+ { dSP; sv_setuv(svarg, beg); GvSV(PL_defgv) = svarg; PUSHMARK(SP); call_sv((SV*)cv, G_VOID|G_DISCARD); }
+ beg += 1 + (beg > 2);
+ }
+ while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
+ START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base )
+ { dSP; sv_setuv(svarg, seg_base + p); GvSV(PL_defgv) = svarg; PUSHMARK(SP); call_sv((SV*)cv, G_VOID|G_DISCARD); }
+ END_DO_FOR_EACH_SIEVE_PRIME
+ }
}
+ end_segment_primes(ctx);
SvREFCNT_dec(svarg);
#endif
XSRETURN_UNDEF;
diff --git a/sieve.c b/sieve.c
index a8dae1a..5e94d5d 100644
--- a/sieve.c
+++ b/sieve.c
@@ -293,3 +293,100 @@ int sieve_segment(unsigned char* mem, UV startd, UV endd)
release_prime_cache(sieve);
return 1;
}
+
+/**************************************************************************/
+
+typedef struct {
+ UV lod;
+ UV hid;
+ UV low;
+ UV high;
+ UV endp;
+ UV segment_size;
+ unsigned char* segment;
+ unsigned char* base;
+} segment_context_t;
+
+/*
+ * unsigned char* segment;
+ * UV seg_base, seg_low, seg_high;
+ * void* ctx = start_segment_primes(low, high, &segment);
+ * while (beg < 7) {
+ * beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
+ * .... with beg ....
+ * beg += 1 + (beg > 2);
+ * }
+ * while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
+ * START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base )
+ * .... with seg_base + p ....
+ * END_DO_FOR_EACH_SIEVE_PRIME
+ * }
+ * end_segment_primes(ctx);
+ */
+
+void* start_segment_primes(UV low, UV high, unsigned char** segmentmem)
+{
+ segment_context_t* ctx;
+
+ MPUassert( high >= low, "start_segment_primes bad arguments");
+ New(0, ctx, 1, segment_context_t);
+ ctx->low = low;
+ ctx->high = high;
+ ctx->lod = low / 30;
+ ctx->hid = high / 30;
+ ctx->endp = (ctx->hid >= (UV_MAX/30)) ? UV_MAX-2 : 30*ctx->hid+29;
+
+ ctx->segment = get_prime_segment( &(ctx->segment_size) );
+ if (ctx->segment == 0)
+ croak("start_segment_primes: Could not get segment");
+ *segmentmem = ctx->segment;
+
+ ctx->base = sieve_erat30( isqrt(ctx->endp)+1 );
+ if (ctx->base == 0)
+ croak("start_segment_primes: Could not get base");
+
+ return (void*) ctx;
+}
+
+int next_segment_primes(void* vctx, UV* base, UV* low, UV* high)
+{
+ UV seghigh_d, range_d;
+ segment_context_t* ctx = (segment_context_t*) vctx;
+
+ if (ctx->lod > ctx->hid) return 0;
+
+ seghigh_d = ((ctx->hid - ctx->lod) < ctx->segment_size)
+ ? ctx->hid
+ : (ctx->lod + ctx->segment_size - 1);
+ range_d = seghigh_d - ctx->lod + 1;
+ *low = ctx->low;
+ *high = (seghigh_d == ctx->hid) ? ctx->high : (seghigh_d*30 + 29);
+ *base = ctx->lod * 30;
+
+ MPUassert( seghigh_d >= ctx->lod, "next_segment_primes: highd < lowd");
+ MPUassert( range_d <= ctx->segment_size, "next_segment_primes: range > segment size");
+
+ if (sieve_segment(ctx->segment, ctx->lod, seghigh_d) == 0) {
+ croak("Could not segment sieve from %"UVuf" to %"UVuf, *base+1, *high);
+ }
+
+ ctx->lod += range_d;
+ ctx->low = *high + 2;
+
+ return 1;
+}
+
+void end_segment_primes(void* vctx)
+{
+ segment_context_t* ctx = (segment_context_t*) vctx;
+ MPUassert(ctx != 0, "end_segment_primes given a null pointer");
+ if (ctx->segment != 0) {
+ release_prime_segment(ctx->segment);
+ ctx->segment = 0;
+ }
+ if (ctx->base != 0) {
+ Safefree(ctx->base);
+ ctx->base = 0;
+ }
+ Safefree(ctx);
+}
diff --git a/sieve.h b/sieve.h
index 4dddde0..1bbc0b9 100644
--- a/sieve.h
+++ b/sieve.h
@@ -6,6 +6,9 @@
extern unsigned char* sieve_erat30(UV end);
extern int sieve_segment(unsigned char* mem, UV startd, UV endd);
+extern void* start_segment_primes(UV low, UV high, unsigned char** segmentmem);
+extern int next_segment_primes(void* vctx, UV* base, UV* low, UV* high);
+extern void end_segment_primes(void* vctx);
static const UV wheel30[] = {1, 7, 11, 13, 17, 19, 23, 29};
diff --git a/t/32-iterators.t b/t/32-iterators.t
index 202ad66..bfc285b 100644
--- a/t/32-iterators.t
+++ b/t/32-iterators.t
@@ -8,7 +8,8 @@ use Math::Prime::Util qw/primes forprimes prime_iterator/;
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
plan tests => 8 + 5
- + 3 + 7;
+ + 3 + 7
+ + 2;
ok(!eval { forprimes { 1 } undef; }, "forprimes undef");
ok(!eval { forprimes { 1 } 2, undef; }, "forprimes 2,undef");
@@ -62,3 +63,28 @@ ok(!eval { prime_iterator(4.5); }, "iterator 4.5");
{my $it = prime_iterator(31398);
is_deeply( [map { $it->() } 1..3], [31469,31477,31481], "iterator 3 primes starting at 31398" );
}
+
+# For fun, nest them.
+{
+ my @t;
+ forprimes {
+ forprimes {
+ forprimes { push @t, $_ } $_,$_+10;
+ } 10*$_,10*$_+10;
+ } 10;
+ is_deeply( [@t], [qw/23 29 31 29 31 37 31 37 41 37 41 43 47 53 59 61 59 61 67 71 73 79 73 79 83 79 83 89/], "triple nested forprimes" );
+}
+{
+ my @t;
+ my $ita = prime_iterator();
+ while ((my $a = $ita->()) <= 10) {
+ my $itb = prime_iterator(10*$a);
+ while ((my $b = $itb->()) <= 10*$a+10) {
+ my $itc = prime_iterator($b);
+ while ((my $c = $itc->()) <= $b+10) {
+ push @t, $c;
+ }
+ }
+ }
+ is_deeply( [@t], [qw/23 29 31 29 31 37 31 37 41 37 41 43 47 53 59 61 59 61 67 71 73 79 73 79 83 79 83 89/], "triple nested iterator" );
+}
diff --git a/util.c b/util.c
index 3507eaa..f0742f4 100644
--- a/util.c
+++ b/util.c
@@ -916,6 +916,8 @@ IV _XS_mertens(UV n) {
#endif
}
+/* TODO: the two Chebyshev functions ought to use segmented sieves */
+
double _XS_chebyshev_theta(UV n)
{
KAHAN_INIT(sum);
--
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