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 f7169a11b7a4d57c7c08bf4434f912f9bbadc39b
Author: Dana Jacobsen <d...@acm.org>
Date: Sat Jun 9 05:54:48 2012 -0600
Overflow for nth_prime, segment prime_count
---
Changes | 3 ++
TODO | 4 +-
XS.xs | 18 ++++++--
t/13-primecount.t | 68 ++++++++++++++++++---------
t/14-nthprime.t | 16 ++++++-
util.c | 135 +++++++++++++++++++++++++++++++++++++++++++++---------
util.h | 1 +
7 files changed, 197 insertions(+), 48 deletions(-)
diff --git a/Changes b/Changes
index 980499a..eaf0d16 100644
--- a/Changes
+++ b/Changes
@@ -2,6 +2,9 @@ Revision history for Perl extension Math::Prime::Util.
0.05 8 June 2012
- Use assembler for mulmod if 64-bit and gcc and x86
+ - nth_prime routines should now all croak on overflow in the same way.
+ - Segmented prime_count, things like this are efficient:
+ say prime_count( 10**16, 10**16 + 2**20 )
0.04 7 June 2012
- Didn't do tests on 32-bit machine before release. Test suite caught
diff --git a/TODO b/TODO
index 8644218..36bdc9e 100644
--- a/TODO
+++ b/TODO
@@ -21,6 +21,6 @@
- speed up random_prime for large numbers
-- Add other bases to pseudoprime test
-
- test x64 asm if Perl is built for 32-bit
+
+- Clean up prime_count ranges. Write tests.
diff --git a/XS.xs b/XS.xs
index 0188258..a73f3c2 100644
--- a/XS.xs
+++ b/XS.xs
@@ -21,13 +21,19 @@ void
prime_memfree()
UV
-prime_count(IN UV n)
+prime_count(IN UV low, IN UV high = 0)
CODE:
+ if (high == 0) { /* Without a Perl layer in front of this, we'll have */
+ high = low; /* the pathological case of a-0 turning into 0-a. */
+ low = 0;
+ }
if (GIMME_V == G_VOID) {
- prime_precalc(n);
+ prime_precalc(high);
RETVAL = 0;
+ } else if (low == 0) {
+ RETVAL = prime_count(high);
} else {
- RETVAL = prime_count(n);
+ RETVAL = prime_count_seg(low, high);
}
OUTPUT:
RETVAL
@@ -140,6 +146,12 @@ segment_primes(IN UV low, IN UV high, IN UV segment_size =
65536UL)
/* To protect vs. overflow, work entirely with d. */
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
diff --git a/t/13-primecount.t b/t/13-primecount.t
index a0711b7..789b259 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -8,30 +8,45 @@ 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 => 10*3 + ($extra ? 10 : 7) + ($use64 ? 2*9 : 0);
+plan tests => 14*3 + ($extra ? 14 : 8) + ($use64 ? 2*18 : 0);
+# Powers of 2: http://oeis.org/A007053/b007053.txt
+# Powers of 10: http://oeis.org/A006880/b006880.txt
my %pivals32 = (
- 1 => 0,
- 10 => 4,
- 100 => 25,
- 1000 => 168,
- 10000 => 1229,
- 100000 => 9592,
- 1000000 => 78498,
- 10000000 => 664579,
- 100000000 => 5761455,
- 1000000000 => 50847534,
+ 1 => 0,
+ 10 => 4,
+ 100 => 25,
+ 1000 => 168,
+ 10000 => 1229,
+ 100000 => 9592,
+ 1000000 => 78498,
+ 10000000 => 664579,
+ 100000000 => 5761455,
+ 1000000000 => 50847534,
+ 65535 => 6542,
+ 16777215 => 1077871,
+ 2147483647 => 105097565,
+ 4294967295 => 203280221,
);
my %pivals64 = (
- 10000000000 => 455052511,
- 100000000000 => 4118054813,
- 1000000000000 => 37607912018,
- 10000000000000 => 346065536839,
- 100000000000000 => 3204941750802,
- 1000000000000000 => 29844570422669,
- 10000000000000000 => 279238341033925,
- 100000000000000000 => 2623557157654233,
-1000000000000000000 => 24739954287740860,
+ 10000000000 => 455052511,
+ 100000000000 => 4118054813,
+ 1000000000000 => 37607912018,
+ 10000000000000 => 346065536839,
+ 100000000000000 => 3204941750802,
+ 1000000000000000 => 29844570422669,
+ 10000000000000000 => 279238341033925,
+ 100000000000000000 => 2623557157654233,
+ 1000000000000000000 => 24739954287740860,
+10000000000000000000 => 234057667276344607,
+ 68719476735 => 2874398515,
+ 1099511627775 => 41203088796,
+ 17592186044415 => 597116381732,
+ 281474976710655 => 8731188863470,
+ 4503599627370495 => 128625503610475,
+ 72057594037927935 => 1906879381028850,
+ 1152921504606846975 => 28423094496953330,
+18446744073709551615 => 425656284035217743,
);
while (my($n, $pin) = each (%pivals32)) {
cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
@@ -40,7 +55,7 @@ while (my($n, $pin) = each (%pivals32)) {
is( prime_count($n), $pin, "Pi($n) = $pin" );
}
my $approx_range = abs($pin - prime_count_approx($n));
- my $range_limit = 1100;
+ my $range_limit = ($n <= 1000000000) ? 1100 : 70000;
cmp_ok( $approx_range, '<=', $range_limit, "prime_count_approx($n) within
$range_limit");
}
if ($use64) {
@@ -50,3 +65,14 @@ if ($use64) {
}
}
+
+# TODO: intervals. From primesieve:
+# 155428406, // prime count 2^32 interval starting at 10^12
+# 143482916, // prime count 2^32 interval starting at 10^13
+# 133235063, // prime count 2^32 interval starting at 10^14
+# 124350420, // prime count 2^32 interval starting at 10^15
+# 116578809, // prime count 2^32 interval starting at 10^16
+# 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 1164387..8f35691 100644
--- a/t/14-nthprime.t
+++ b/t/14-nthprime.t
@@ -8,7 +8,7 @@ use Math::Prime::Util qw/nth_prime nth_prime_lower
nth_prime_upper nth_prime_app
my $use64 = Math::Prime::Util::_maxbits > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
-plan tests => 7*2 + 9*3 + ($extra ? 9 : 7) + ($use64 ? 9*3 : 0);
+plan tests => 7*2 + 9*3 + 7 + ($extra ? 9 : 7) + ($use64 ? 9*3 : 0);
my %pivals32 = (
1 => 0,
@@ -26,6 +26,8 @@ while (my($n, $pin) = each (%pivals32)) {
cmp_ok( nth_prime($next), '>=', $n, "nth_prime($next) >= $n");
}
+# Powers of 10: http://oeis.org/A006988/b006988.txt
+# Powers of 2: http://oeis.org/A033844/b033844.txt
my %nthprimes32 = (
1 => 2,
10 => 29,
@@ -73,3 +75,15 @@ if ($use64) {
}
}
+my $maxindex = $use64 ? 425656284035217743 : 203280221;
+my $maxprime = $use64 ? 18446744073709551557 : 4294967291;
+cmp_ok( nth_prime_lower($maxindex), '<=', $maxprime,
"nth_prime_lower(maxindex) <= maxprime");
+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");
diff --git a/util.c b/util.c
index f969af5..c924322 100644
--- a/util.c
+++ b/util.c
@@ -284,7 +284,7 @@ static const double F1 = 1.0;
UV prime_count_lower(UV x)
{
double fx, flogx;
- double a = 1.80;
+ double a = 1.80; /* Dusart 1999, page 14 */
if (x < NPRIME_COUNT_SMALL)
return prime_count_small[x];
@@ -314,7 +314,7 @@ UV prime_count_lower(UV x)
UV prime_count_upper(UV x)
{
double fx, flogx;
- double a = 2.51;
+ double a = 2.51; /* Dusart 1999, page 14 */
if (x < NPRIME_COUNT_SMALL)
return prime_count_small[x];
@@ -471,42 +471,49 @@ static const unsigned short primes_small[] =
409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499};
#define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))
-/* The nth prime will be more than this number */
+/* The nth prime will be greater than or equal to this number */
UV nth_prime_lower(UV n)
{
double fn = (double) n;
double flogn, flog2n, lower;
if (n < NPRIMES_SMALL)
- return (n==0) ? 0 : primes_small[n]-1;
+ return (n==0) ? 0 : primes_small[n];
flogn = log(n);
- flog2n = log(flogn);
+ flog2n = log(flogn); /* Note distinction between log_2(n) and log^2(n) */
- /* Dusart 1999, for all n >= 2 */
+ /* Dusart 1999 page 14, for all n >= 2 */
lower = fn * (flogn + flog2n - 1.0 + ((flog2n-2.25)/flogn));
- if (lower > (double)UV_MAX)
- return 0;
+ /* Watch out for overflow */
+ if (lower >= (double)UV_MAX) {
+#if BITS_PER_WORD == 32
+ if (n <= UVCONST(203280221)) return UVCONST(4294967291);
+#else
+ if (n <= UVCONST(425656284035217743)) return UVCONST(18446744073709551557);
+#endif
+ croak("nth_prime_lower(%"UVuf") overflow", n);
+ }
return (UV) lower;
}
-/* The nth prime will be less than this number */
+/* The nth prime will be less or equal to this number */
UV nth_prime_upper(UV n)
{
double fn = (double) n;
double flogn, flog2n, upper;
if (n < NPRIMES_SMALL)
- return primes_small[n]+1;
+ return primes_small[n];
flogn = log(n);
- flog2n = log(flogn);
+ flog2n = log(flogn); /* Note distinction between log_2(n) and log^2(n) */
if (n >= 39017)
- upper = fn * ( flogn + flog2n - 0.9484 ); /* Dusart 1999 */
+ upper = fn * ( flogn + flog2n - 0.9484 ); /* Dusart 1999 page 14*/
else if (n >= 27076)
upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.80)/flogn)); /*Dusart
1999*/
else if (n >= 7022)
@@ -514,9 +521,15 @@ UV nth_prime_upper(UV n)
else
upper = fn * ( flogn + flog2n );
- /* Special case to not overflow any 32-bit */
- if (upper > (double)UV_MAX)
- return (n <= UVCONST(203280221)) ? UVCONST(4294967292) : 0;
+ /* Watch out for overflow */
+ if (upper >= (double)UV_MAX) {
+#if BITS_PER_WORD == 32
+ if (n <= UVCONST(203280221)) return UVCONST(4294967291);
+#else
+ if (n <= UVCONST(425656284035217743)) return UVCONST(18446744073709551557);
+#endif
+ croak("nth_prime_upper(%"UVuf") overflow", n);
+ }
return (UV) ceil(upper);
}
@@ -542,6 +555,8 @@ UV nth_prime_approx(UV n)
* m=1 + ((flog2n - 2)/flogn) );
* m=2 - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
* + O((flog2n/flogn)^3)
+ *
+ * Shown in Dusart 1999 page 12, as well as other sources
*/
approx = fn * ( flogn + flog2n - 1
+ ((flog2n - 2)/flogn)
@@ -560,8 +575,24 @@ UV nth_prime_approx(UV n)
else if (n < 12000) approx += 3.0 * order;
else if (n <150000) approx += 2.1 * order;
- if (approx > (double)UV_MAX)
- return 0;
+ /* For all three analytical functions, it is possible that for a given valid
+ * input, we will not be able to return an output that fits in the UV type.
+ * For example, if they ask for the 203280222nd prime, we should return
+ * 4294967311. But in 32-bit, that overflows. What we do is calculate our
+ * double precision value. If that would overflow, then we look at the input
+ * and if it is <= the index of the last representable prime, then we return
+ * the last representable prime. Otherwise, we croak an overflow message.
+ * This should maintain the invariant:
+ * nth_prime_lower(n) <= nth_prime(n) <= nth_prime_upper(n)
+ */
+ if (approx >= (double)UV_MAX) {
+#if BITS_PER_WORD == 32
+ if (n <= UVCONST(203280221)) return UVCONST(4294967291);
+#else
+ if (n <= UVCONST(425656284035217743)) return UVCONST(18446744073709551557);
+#endif
+ croak("nth_prime_approx(%"UVuf") overflow", n);
+ }
return (UV) (approx + 0.5);
}
@@ -622,10 +653,7 @@ UV nth_prime(UV n)
/* Determine a bound on the nth prime. We know it comes before this. */
upper_limit = nth_prime_upper(n);
- if (upper_limit == 0) {
- croak("nth_prime(%"UVuf") would overflow", n);
- return 0;
- }
+ MPUassert(upper_limit > 0, "nth_prime got an upper limit of 0");
/* Get the primary cache, and ensure it is at least this large. If the
* input is small enough, get a sieve covering the range. Otherwise, we'll
@@ -669,3 +697,68 @@ UV nth_prime(UV n)
MPUassert(count == target, "nth_prime got incorrect count");
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;
+}
diff --git a/util.h b/util.h
index b14c175..6a09d98 100644
--- a/util.h
+++ b/util.h
@@ -13,6 +13,7 @@ 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 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
_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
