[libmath-prime-util-perl] 09/16: Rewrite prime_count for segments
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:25 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.06
in repository libmath-prime-util-perl.
commit cd97a8605a0a36236383719ab7e0da19342ae1e2
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Jun 11 03:43:08 2012 -0600
Rewrite prime_count for segments
---
TODO | 6 +-
XS.xs | 4 +-
sieve.c | 6 +-
t/13-primecount.t | 50 +++++++-
t/14-nthprime.t | 32 +++--
t/17-pseudoprime.t | 7 +-
t/18-functions.t | 4 +-
util.c | 355 ++++++++++++++++++++++++++++-------------------------
util.h | 3 +-
9 files changed, 267 insertions(+), 200 deletions(-)
diff --git a/TODO b/TODO
index 36bdc9e..2661193 100644
--- a/TODO
+++ b/TODO
@@ -9,18 +9,16 @@
- Add test to check maxbits in compiled library vs. Perl
-- Li(n)
-
- Pure perl implementations
- input validation (in XS, or do we need to make Perl wrappers for everything?)
- Faster SQUFOF
-- nth_prime bounds are < in the code, but <= in the documentation.
-
- speed up random_prime for large numbers
- test x64 asm if Perl is built for 32-bit
- Clean up prime_count ranges. Write tests.
+
+- better prime count upper/lower bounds
diff --git a/XS.xs b/XS.xs
index e12bbe5..bf9a5b7 100644
--- a/XS.xs
+++ b/XS.xs
@@ -30,10 +30,8 @@ prime_count(IN UV low, IN UV high = 0)
if (GIMME_V == G_VOID) {
prime_precalc(high);
RETVAL = 0;
- } else if (low == 0) {
- RETVAL = prime_count(high);
} else {
- RETVAL = prime_count_seg(low, high);
+ RETVAL = prime_count(low, high);
}
OUTPUT:
RETVAL
diff --git a/sieve.c b/sieve.c
index 6ac1af6..0636f4b 100644
--- a/sieve.c
+++ b/sieve.c
@@ -32,8 +32,10 @@ UV get_prime_cache(UV n, const unsigned char** sieve)
prime_cache_size = 0;
/* Sieve a bit more than asked, to mitigate thrashing */
- if (n < (UV_MAX-3840))
- n += 3840;
+ if (n >= (UV_MAX-3840))
+ n = UV_MAX;
+ else
+ n = ((n + 3840)/30)*30;
/* TODO: testing near 2^32-1 */
prime_cache_sieve = sieve_erat30(n);
diff --git a/t/13-primecount.t b/t/13-primecount.t
index 789b259..067e786 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/prime_count prime_count_lower prime_count_upper prime_c
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-plan tests => 14*3 + ($extra ? 14 : 8) + ($use64 ? 2*18 : 0);
+plan tests => 14*3 + 8 + 6*$extra + 2*18*$use64 + 12 + 5*$use64;
# Powers of 2: http://oeis.org/A007053/b007053.txt
# Powers of 10: http://oeis.org/A006880/b006880.txt
@@ -65,7 +65,52 @@ if ($use64) {
}
}
+# ./primesieve 1e10 -o2**32 -c1
+# ./primesieve 24689 7973249 -c1
+my %intervals = (
+ "868396 to 9478505" => 563275,
+ "1118105 to 9961674" => 575195,
+ "24689 to 7973249" => 535368,
+ "1e10 +2**16" => 2821,
+ "17 to 13" => 0,
+ "3 to 17" => 6,
+ "4 to 17" => 5,
+ "4 to 16" => 4,
+ "191912783 +248" => 2,
+ "191912784 +247" => 1,
+ "191912783 +247" => 1,
+ "191912784 +246" => 0,
+ "1e14 +2**16" => 1973,
+ "127976334671 +468" => 2,
+ "127976334672 +467" => 1,
+ "127976334671 +467" => 1,
+ "127976334672 +466" => 0,
+);
+
+while (my($range, $expect) = each (%intervals)) {
+ my($low,$high);
+ if ($range =~ /(\S+)\s+to\s+(\S+)/) {
+ $low = fixnum($1);
+ $high = fixnum($2);
+ } elsif ($range =~ /(\S+)\s*\+\s*(\S+)/) {
+ $low = fixnum($1);
+ $high = $low + fixnum($2);
+ } else {
+ die "Can't parse test data";
+ }
+ next if $high > ~0;
+ is( prime_count($low,$high), $expect, "prime_count($range) = $expect");
+}
+
+sub fixnum {
+ my $nstr = shift;
+ $nstr =~ s/^(\d+)e(\d+)$/$1*(10**$2)/e;
+ $nstr =~ s/^(\d+)\*\*(\d+)$/$1**$2/e;
+ die "Unknown string in test" unless $nstr =~ /^\d+$/;
+ $nstr;
+}
+
# TODO: intervals. From primesieve:
# 155428406, // prime count 2^32 interval starting at 10^12
# 143482916, // prime count 2^32 interval starting at 10^13
@@ -75,4 +120,5 @@ if ($use64) {
# 109726486, // prime count 2^32 interval starting at 10^17
# 103626726, // prime count 2^32 interval starting at 10^18
# 98169972}; // prime count 2^32 interval starting at 10^19
-# Note: ./primesieve 1e10 -o2**32 -c1
+
+
diff --git a/t/14-nthprime.t b/t/14-nthprime.t
index 8f35691..e777ca9 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -3,12 +3,16 @@ use strict;
use warnings;
use Test::More;
-use Math::Prime::Util qw/nth_prime nth_prime_lower nth_prime_upper nth_prime_approx/;
+use Math::Prime::Util qw/primes nth_prime nth_prime_lower nth_prime_upper nth_prime_approx/;
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
+my $broken64 = (18446744073709550592 == ~0);
-plan tests => 7*2 + 9*3 + 7 + ($extra ? 9 : 7) + ($use64 ? 9*3 : 0);
+my $nsmallprimes = 1000;
+my $nth_small_prime = 7919; # nth_prime(1000)
+
+plan tests => 7*2 + $nsmallprimes+1 + 9*3 + 7 + ($extra ? 9 : 7) + ($use64 ? 9*3 : 0);
my %pivals32 = (
1 => 0,
@@ -26,6 +30,13 @@ while (my($n, $pin) = each (%pivals32)) {
cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
}
+my @small_primes = (0);
+push @small_primes, @{primes($nth_small_prime)};
+foreach my $n (0 .. $nsmallprimes) {
+ is(nth_prime($n), $small_primes[$n], "The ${n}th prime is $small_primes[$n]");
+}
+
+
# Powers of 10: http://oeis.org/A006988/b006988.txt
# Powers of 2: http://oeis.org/A033844/b033844.txt
my %nthprimes32 = (
@@ -81,9 +92,14 @@ cmp_ok( nth_prime_lower($maxindex), '<=', $maxprime, "nth_prime_lower(maxindex)
cmp_ok( nth_prime_upper($maxindex), '>=', $maxprime, "nth_prime_upper(maxindex) >= maxprime");
cmp_ok( nth_prime_approx($maxindex), '==', $maxprime, "nth_prime_approx(maxindex) == maxprime");
cmp_ok( nth_prime_lower($maxindex+1), '>=', nth_prime_lower($maxindex), "nth_prime_lower(maxindex+1) >= nth_prime_lower(maxindex)");
-eval { nth_prime_upper($maxindex+1); };
-like($@, qr/overflow/, "nth_prime_upper(maxindex+1) overflows");
-eval { nth_prime_approx($maxindex+1); };
-like($@, qr/overflow/, "nth_prime_approx(maxindex+1) overflows");
-eval { nth_prime($maxindex+1); };
-like($@, qr/overflow/, "nth_prime(maxindex+1) overflows");
+
+my $add = ($broken64) ? 100 : 1;
+
+eval { nth_prime_upper($maxindex+$add); };
+like($@, qr/overflow/, "nth_prime_upper(maxindex+$add) overflows");
+
+eval { nth_prime_approx($maxindex+$add); };
+like($@, qr/overflow/, "nth_prime_approx(maxindex+$add) overflows");
+
+eval { nth_prime($maxindex+$add); };
+like($@, qr/overflow/, "nth_prime(maxindex+$add) overflows");
diff --git a/t/17-pseudoprime.t b/t/17-pseudoprime.t
index bc2de62..1c7b18c 100644
--- a/t/17-pseudoprime.t
+++ b/t/17-pseudoprime.t
@@ -8,10 +8,7 @@ use Math::Prime::Util qw/is_prime miller_rabin/;
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-#plan tests => 6 + 19 + 3573 + (5 + 29 + 22 + 23 + 16) + 15 + 27
-# + ($use64 ? 5+1 : 0)
-# + ($extra ? 6 : 0)
-# + (($extra && $use64) ? 19 : 0);
+plan tests => 3 + 295 + 4 + 4*$use64 + 1 + 1*$extra + 161;
ok(!eval { miller_rabin(2047); }, "MR with no base fails");
ok(!eval { miller_rabin(2047,0); }, "MR base 0 fails");
@@ -105,5 +102,3 @@ foreach my $base (@ebases) {
# 1005905886, 1340600841, 553174392, 3046413974,
# 203659041, 3613982119,
# 9780504, 1795265022
-
-done_testing();
diff --git a/t/18-functions.t b/t/18-functions.t
index 1ed634a..4e5bda7 100644
--- a/t/18-functions.t
+++ b/t/18-functions.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/prime_count ExponentialIntegral LogarithmicIntegral Rie
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-#plan tests => 14*3 + ($extra ? 14 : 8) + ($use64 ? 2*18 : 0);
+plan tests => 4 + 1 + 12 + 11 + 9;
eval { ExponentialIntegral(0); };
like($@, qr/invalid/i, "Ei(0) is invalid");
@@ -77,8 +77,6 @@ while (my($n, $rin) = each (%rvals)) {
}
-done_testing();
-
sub cmp_closeto {
my $got = shift;
my $expect = shift;
diff --git a/util.c b/util.c
index c469f05..ab65a0c 100644
--- a/util.c
+++ b/util.c
@@ -3,6 +3,11 @@
#include <string.h>
#include <math.h>
+/* It took until C99 to get this macro */
+#ifndef INFINITY
+#define INFINITY (DBL_MAX + DBL_MAX)
+#endif
+
#include "util.h"
#include "sieve.h"
#include "factor.h"
@@ -237,6 +242,109 @@ UV prev_prime(UV n)
+/* Given a sieve of size nbytes, walk it counting zeros (primes) until:
+ *
+ * (1) we counted them all: return the count, which will be less than maxcount.
+ *
+ * (2) we hit maxcount: set position to the index of the maxcount'th prime
+ * and return count (which will be equal to maxcount).
+ */
+static UV count_segment_maxcount(const unsigned char* sieve, UV nbytes, UV maxcount, UV* pos)
+{
+ UV count = 0;
+ UV byte = 0;
+
+ MPUassert(sieve != 0, "count_segment_maxcount incorrect args");
+ MPUassert(pos != 0, "count_segment_maxcount incorrect args");
+ *pos = 0;
+ if ( (nbytes == 0) || (maxcount == 0) )
+ return 0;
+
+ while ( (byte < nbytes) && (count < maxcount) )
+ count += byte_zeros[sieve[byte++]];
+
+ if (count >= maxcount) { /* One too far -- back up */
+ count -= byte_zeros[sieve[--byte]];
+ }
+
+ MPUassert(count < maxcount, "count_segment_maxcount wrong count");
+
+ if (byte == nbytes)
+ return count;
+
+ /* The result is somewhere in the next byte */
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, byte*30+1, nbytes*30-1)
+ if (++count == maxcount) { *pos = p; return count; }
+ END_DO_FOR_EACH_SIEVE_PRIME;
+
+ MPUassert(0, "count_segment_maxcount failure");
+ return 0;
+}
+
+/* Given a sieve of size nbytes, counting zeros (primes) but excluding the
+ * areas outside lowp and highp.
+ */
+static UV count_segment_ranged(const unsigned char* sieve, UV nbytes, UV lowp, UV highp)
+{
+ UV count = 0;
+
+ UV lo_d = lowp/30;
+ UV lo_m = lowp - lo_d*30;
+ UV hi_d = highp/30;
+ UV hi_m = highp - hi_d*30;
+
+ MPUassert( sieve != 0, "count_segment_maxcount incorrect args");
+
+ if (hi_d >= nbytes) {
+ hi_d = nbytes-1;
+ highp = hi_d*30+29;
+ }
+
+ if ( (nbytes == 0) || (highp < lowp) )
+ return 0;
+
+#if 0
+ /* Dead simple way */
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, lowp, highp)
+ count++;
+ END_DO_FOR_EACH_SIEVE_PRIME;
+ return count;
+#endif
+
+ /* Count first fragment */
+ if (lo_m > 1) {
+ UV upper = (highp <= (lo_d*30+29)) ? highp : (lo_d*30+29);
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, lowp, upper)
+ count++;
+ END_DO_FOR_EACH_SIEVE_PRIME;
+ lowp = upper+2;
+ lo_d = lowp/30;
+ lo_m = lowp - lo_d*30;
+ }
+ if (highp < lowp)
+ return count;
+
+ /* Count bytes in the middle */
+ {
+ UV count_bytes = hi_d - lo_d + (hi_m == 29);
+ if (count_bytes > 0) {
+ count += count_zero_bits(sieve+lo_d, count_bytes);
+ lowp += 30*count_bytes;
+ lo_d = lowp/30;
+ lo_m = lowp - lo_d*30;
+ }
+ }
+ if (highp < lowp)
+ return count;
+
+ /* Count last fragment */
+ START_DO_FOR_EACH_SIEVE_PRIME(sieve, lowp, highp)
+ count++;
+ END_DO_FOR_EACH_SIEVE_PRIME;
+
+ return count;
+}
+
/*
* The pi(x) prime count functions. prime_count(x) gives an exact number,
@@ -295,17 +403,23 @@ UV prime_count_lower(UV x)
if (x < 599)
return (UV) (fx / (flogx-0.7));
- if (x < 2700) { a = 0.30; }
- else if (x < 5500) { a = 0.90; }
- else if (x < 19400) { a = 1.30; }
- else if (x < 32299) { a = 1.60; }
- else if (x < 176000) { a = 1.80; }
- else if (x < 315000) { a = 2.10; }
- else if (x < 1100000) { a = 2.20; }
- else if (x < 4500000) { a = 2.31; }
- else if (x <233000000) { a = 2.36; }
- else if (x <240000000) { a = 2.32; }
- else if (x <UVCONST(0xFFFFFFFF)) { a = 2.32; }
+ if (x < 2700) { a = 0.30; }
+ else if (x < 5500) { a = 0.90; }
+ else if (x < 19400) { a = 1.30; }
+ else if (x < 32299) { a = 1.60; }
+ else if (x < 176000) { a = 1.80; }
+ else if (x < 315000) { a = 2.10; }
+ else if (x < 1100000) { a = 2.20; }
+ else if (x < 4500000) { a = 2.31; }
+ else if (x <233000000) { a = 2.36; }
+ else if (x <240000000) { a = 2.32; }
+#if BITS_PER_WORD == 32
+ else a = 2.32;
+#else
+ else if (x < UVCONST( 5433800000)) { a = 2.32; }
+ else if (x < UVCONST(10000000000)) { a = 2.15; }
+ else a = 1.80;
+#endif
return (UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) );
}
@@ -351,7 +465,12 @@ UV prime_count_upper(UV x)
else if (x <400000000) { a = 2.375; }
else if (x <510000000) { a = 2.370; }
else if (x <682000000) { a = 2.367; }
- else if (x <UVCONST(0xFFFFFFFF)) { a = 2.362; }
+#if BITS_PER_WORD == 32
+ else a = 2.362;
+#else
+ else if (x < UVCONST(10000000000)) { a = 2.362; }
+ else a = 2.51;
+#endif
/*
* An alternate idea:
@@ -404,75 +523,73 @@ UV prime_count_approx(UV x)
}
-UV prime_count(UV n)
+UV prime_count(UV low, UV high)
{
- const unsigned char* sieve;
- static UV last_bytes = 0;
- static UV last_count = 3;
- UV s, bytes;
- UV count = 3;
-
- if (n < NPRIME_COUNT_SMALL)
- return prime_count_small[n];
-
- bytes = n/30;
- s = 0;
-
- if (n <= get_prime_cache(0, &sieve)) {
-
- /* We have enough primes -- just count them. */
-
- /* Start from last word position if we can. This is a big speedup when
- * calling prime_count many times with successively larger numbers. */
- if (bytes >= last_bytes) {
- s = last_bytes;
- count = last_count;
- }
+ const unsigned char* cache_sieve;
+ unsigned char* segment;
+ UV segbase, segment_size, low_d, high_d;
+ UV count = 0;
- count += count_zero_bits(sieve+s, bytes-s);
+ if ( (low <= 2) && (high < NPRIME_COUNT_SMALL) )
+ return prime_count_small[high];
- last_bytes = bytes;
- last_count = count;
+ if ((low <= 2) && (high >= 2)) count++;
+ if ((low <= 3) && (high >= 3)) count++;
+ if ((low <= 5) && (high >= 5)) count++;
+ if (low < 7) low = 7;
- START_DO_FOR_EACH_SIEVE_PRIME(sieve, 30*bytes, n)
- count++;
- END_DO_FOR_EACH_SIEVE_PRIME;
+ if (low > high) return count;
- } else {
+ low_d = low/30;
+ high_d = high/30;
+ //printf(" our range %lu - %lu lies in bytes %lu - %lu\n", low, high, low_d, high_d);
- /* We don't have enough primes. Repeatedly segment sieve */
- UV const segment_size = SEGMENT_CHUNK_SIZE;
- unsigned char* segment;
+ /* Count full bytes only -- no fragments from primary cache */
+ segment_size = get_prime_cache(0, &cache_sieve) / 30;
+ if (segment_size < high_d) {
+ /* Expand sieve to sqrt(n) */
+ UV endp = (high_d >= (UV_MAX/30)) ? UV_MAX-2 : 30*high_d+29;
+ segment_size = get_prime_cache( sqrt(endp) + 1 , &cache_sieve) / 30;
+ }
+ //printf(" primary cache is %lu bytes. From %lu to %lu. Size=%lu, count=%lu\n", segment_size, 0, (segment_size-1)*30+29, get_prime_cache_size(), count_zero_bits(cache_sieve, segment_size-1));
- /* Get this over with once */
- prime_precalc( sqrt(n) + 2 );
+ if (1 && low_d <= segment_size) {
+ /* Count all the primes in the primary cache in our range */
+ count += count_segment_ranged(cache_sieve, segment_size, low, high);
+ //printf(" counted in pcache, count (%lu,%lu) is %lu\n", low, high, count);
- segment = get_prime_segment();
- if (segment == 0)
- return 0;
+ if (high_d <= segment_size)
+ return count;
- for (s = 0; s <= bytes; s += segment_size) {
- /* We want to sieve one extra byte, to handle the last fragment */
- UV sieve_bytes = ((bytes-s) >= segment_size) ? segment_size : bytes-s+1;
- UV count_bytes = ((bytes-s) >= segment_size) ? segment_size : bytes-s;
+ low_d = segment_size;
+ }
- /* printf("sieving from %"UVuf" to %"UVuf"\n", 30*s+1, 30*(s+sieve_bytes-1)+29); */
- if (sieve_segment(segment, s, s + sieve_bytes - 1) == 0) {
- croak("Could not segment sieve from %"UVuf" to %"UVuf, 30*s+1, 30*(s+sieve_bytes)+29);
- break;
- }
+ /* More primes needed. Repeatedly segment sieve */
+ segment_size = SEGMENT_CHUNK_SIZE;
+ segment = get_prime_segment();
+ if (segment == 0)
+ return 0;
- if (count_bytes > 0)
- count += count_zero_bits(segment, count_bytes);
+ while (low_d <= high_d)
+ {
+ UV seghigh_d = ((high_d - low_d) < segment_size)
+ ? high_d
+ : (low_d + segment_size-1);
+ UV range_d = seghigh_d - low_d + 1;
+ UV seglow = (low_d*30 < low) ? low : low_d*30;
+ UV seghigh = (seghigh_d == high_d) ? high : (seghigh_d*30+29);
+ //printf(" seg d: %lu - %lu seg p: %lu - %lu\n", low_d, seghigh_d, seglow, seghigh);
+ //printf(" sieve from bytes %lu to %lu\n", low_d, seghigh_d);
+ if (sieve_segment(segment, low_d, seghigh_d) == 0) {
+ croak("Could not segment sieve from %"UVuf" to %"UVuf, segbase+1, 30*(segbase+segment_size)+29);
+ break;
}
- s -= segment_size;
- /*printf("counting fragment from %"UVuf" to %"UVuf"\n", 30*bytes-30*s, n-30*s); */
- START_DO_FOR_EACH_SIEVE_PRIME(segment, 30*bytes - s*30, n - s*30)
- count++;
- END_DO_FOR_EACH_SIEVE_PRIME;
+ //printf(" count from %lu to %lu\n", seglow - low_d*30, seghigh - low_d*30);
+ count += count_segment_ranged(segment, segment_size, seglow - low_d*30, seghigh - low_d*30);
+ low_d += range_d;
}
return count;
@@ -615,46 +732,6 @@ UV nth_prime_approx(UV n)
}
-/*
- * Given a sieve of size nbytes, walk it counting zeros (primes) until:
- *
- * (1) we counted them all: return the count, which will be less than maxcount.
- *
- * (2) we hit maxcount: set position to the index of the maxcount'th prime
- * and return count (which will be equal to maxcount).
- */
-static UV count_segment(const unsigned char* sieve, UV nbytes, UV maxcount, UV* pos)
-{
- UV count = 0;
- UV byte = 0;
-
- MPUassert(sieve != 0, "count_segment incorrect args");
- MPUassert(pos != 0, "count_segment incorrect args");
- *pos = 0;
- if ( (nbytes == 0) || (maxcount == 0) )
- return 0;
-
- while ( (byte < nbytes) && (count < maxcount) )
- count += byte_zeros[sieve[byte++]];
-
- if (count >= maxcount) { /* One too far -- back up */
- count -= byte_zeros[sieve[--byte]];
- }
-
- MPUassert(count < maxcount, "count_segment wrong count");
-
- if (byte == nbytes)
- return count;
-
- /* The result is somewhere in the next byte */
- START_DO_FOR_EACH_SIEVE_PRIME(sieve, byte*30+1, nbytes*30-1)
- if (++count == maxcount) { *pos = p; return count; }
- END_DO_FOR_EACH_SIEVE_PRIME;
-
- MPUassert(0, "count_segment failure");
- return 0;
-}
-
UV nth_prime(UV n)
{
const unsigned char* cache_sieve;
@@ -683,7 +760,7 @@ UV nth_prime(UV n)
segment_size = get_prime_cache(sqrt(upper_limit), &cache_sieve) / 30;
/* Count up everything in the cached sieve. */
- count += count_segment(cache_sieve, segment_size, target, &p);
+ count += count_segment_maxcount(cache_sieve, segment_size, target, &p);
if (count == target)
return p;
@@ -706,7 +783,7 @@ UV nth_prime(UV n)
}
/* Count up everything in this segment */
- count += count_segment(segment, segment_size, target-count, &p);
+ count += count_segment_maxcount(segment, segment_size, target-count, &p);
if (count < target)
segbase += segment_size;
@@ -715,70 +792,8 @@ UV nth_prime(UV n)
return ( (segbase*30) + p );
}
-UV prime_count_seg(UV low, UV high)
-{
- UV const segment_size = SEGMENT_CHUNK_SIZE;
- unsigned char* segment;
- UV low_d, high_d;
- UV count = 0;
-
- if ((low <= 2) && (high >= 2)) count++;
- if ((low <= 3) && (high >= 3)) count++;
- if ((low <= 5) && (high >= 5)) count++;
- if (low < 7) low = 7;
-
- if (low > high) return count;
-
- /*
- * (1) be smart about small low/high values.
- * (2) be generic so we can kill the other function
- * (3) use fast counting.
- */
-
- segment = get_prime_segment();
- if (segment == 0)
- return 0;
-
- low_d = low / 30;
- high_d = high / 30;
-
- { /* Avoid recalculations of this */
- UV endp = (high_d >= (UV_MAX/30)) ? UV_MAX-2 : 30*high_d+29;
- prime_precalc( sqrt(endp) + 1 );
- }
-
- while ( low_d <= high_d ) {
- UV seghigh_d = ((high_d - low_d) < segment_size)
- ? high_d
- : (low_d + segment_size-1);
- UV range_d = seghigh_d - low_d + 1;
- UV seghigh = (seghigh_d == high_d) ? high : (seghigh_d*30+29);
- UV segbase = low_d * 30;
- /* printf(" startd = %"UVuf" endd = %"UVuf"\n", startd, endd); */
-
- MPUassert( seghigh_d >= low_d, "segment_primes highd < lowd");
- MPUassert( range_d <= segment_size, "segment_primes range > segment size");
-
- /* Sieve from startd*30+1 to endd*30+29. */
- if (sieve_segment(segment, low_d, seghigh_d) == 0) {
- croak("Could not segment sieve from %"UVuf" to %"UVuf, segbase+1, seghigh);
- break;
- }
- if (seghigh < high) {
- count += count_zero_bits(segment, range_d);
- } else {
- START_DO_FOR_EACH_SIEVE_PRIME( segment, low - segbase, seghigh - segbase )
- count++;
- END_DO_FOR_EACH_SIEVE_PRIME
- }
- low_d += range_d;
- low = seghigh+2;
- }
-
- return count;
-}
/*
* See:
diff --git a/util.h b/util.h
index cb7ff43..b4eca82 100644
--- a/util.h
+++ b/util.h
@@ -12,8 +12,7 @@ extern UV prev_prime(UV x);
extern UV prime_count_lower(UV x);
extern UV prime_count_upper(UV x);
extern UV prime_count_approx(UV x);
-extern UV prime_count(UV x);
-extern UV prime_count_seg(UV low, UV high);
+extern UV prime_count(UV low, UV high);
extern UV nth_prime_lower(UV n);
extern UV nth_prime_upper(UV x);
--
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