This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.34
in repository libmath-prime-util-perl.
commit 58456992d78b86333cf802636cf44ca63de69a46
Author: Dana Jacobsen <d...@acm.org>
Date: Tue Nov 19 16:16:53 2013 -0800
Switch small primes and lpf arrays to uint32_t
---
Changes | 4 ++
TODO | 3 --
lehmer.c | 76 ++++++++++++++-------------
xt/primecount-many.t | 145 +++++++++++++++++++++++++++++++++++++++++++++++++++
4 files changed, 189 insertions(+), 39 deletions(-)
diff --git a/Changes b/Changes
index 75489cc..796776a 100644
--- a/Changes
+++ b/Changes
@@ -4,6 +4,10 @@ Revision history for Perl module Math::Prime::Util
- Fixed test that was using a 64-bit number on 32-bit machines.
+ - Switch a couple internal arrays from UV to uint32 in prime count.
+ This reduces memory consumption a little with big counts. Total
+ memory use for counts > 10^15 is about 5x less than in version 0.31.
+
0.33 2013-11-18
diff --git a/TODO b/TODO
index 4f2c11d..b336fdd 100644
--- a/TODO
+++ b/TODO
@@ -57,9 +57,6 @@
algorithm). The PP code isn't doing that, which means we're doing lots of
extra primality checks, which aren't cheap in PP.
-- I believe we can make the Lehmer/LMO small primes uint32_t, which will
- give some more memory reduction and a little speed.
-
- use Math::BigInt instead of require, and return bigints as needed.
This is a fundamental change in how some functions operate, though likely
one that most people wouldn't see, and should be less surprise:
diff --git a/lehmer.c b/lehmer.c
index 45f6fe2..0d9501c 100644
--- a/lehmer.c
+++ b/lehmer.c
@@ -142,7 +142,7 @@ static UV isqrt(UV n)
}
/* Callback used for creating an array of primes. */
-static UV* sieve_array = 0;
+static uint32_t* sieve_array = 0;
static UV sieve_k;
static UV sieve_n;
class stop_primesieve : public std::exception { };
@@ -158,10 +158,10 @@ void primesieve_callback(uint64_t pk) {
/* Generate an array of n small primes, where the kth prime is element p[k].
* Remember to free when done. */
#define TINY_PRIME_SIZE 20000
-static UV* tiny_primes = 0;
-static UV* generate_small_primes(UV n)
+static uint32_t* tiny_primes = 0;
+static uint32_t* generate_small_primes(UV n)
{
- UV* primes;
+ uint32_t* primes;
double fn = (double)n;
double flogn = log(fn);
double flog2n = log(flogn);
@@ -170,13 +170,13 @@ static UV* generate_small_primes(UV n)
(n >= 178974) ? (UV) ceil(fn*(flogn+flog2n-1.0+((flog2n-1.95)/flogn))) :
(n >= 18) ? (UV) ceil(fn*(flogn+flog2n-1.0+((flog2n+0.30)/flogn)))
: 59;
- New(0, primes, n+1, UV);
+ New(0, primes, n+1, uint32_t);
if (primes == 0)
croak("Can not allocate small primes\n");
if (n < TINY_PRIME_SIZE) {
if (tiny_primes == 0)
tiny_primes = generate_small_primes(TINY_PRIME_SIZE+1);
- memcpy(primes, tiny_primes, (n+1) * sizeof(UV));
+ memcpy(primes, tiny_primes, (n+1) * sizeof(uint32_t));
return primes;
}
primes[0] = 0;
@@ -205,9 +205,9 @@ static UV* generate_small_primes(UV n)
/* Generate an array of n small primes, where the kth prime is element p[k].
* Remember to free when done. */
-static UV* generate_small_primes(UV n)
+static uint32_t* generate_small_primes(UV n)
{
- UV* primes;
+ uint32_t* primes;
UV i = 0;
double fn = (double)n;
double flogn = log(fn);
@@ -218,7 +218,9 @@ static UV* generate_small_primes(UV n)
(n >= 18) ? (UV) ceil(fn*(flogn+flog2n-1.0+((flog2n+0.30)/flogn)))
: 59;
- New(0, primes, n+1, UV);
+ if (n > 203280221)
+ croak("generate small primes with argument too large: %lu\n", (unsigned
long)n);
+ New(0, primes, n+1, uint32_t);
if (primes == 0)
croak("Can not allocate small primes\n");
primes[0] = 0;
@@ -228,7 +230,7 @@ static UV* generate_small_primes(UV n)
} END_DO_FOR_EACH_PRIME
if (i < n)
croak("Did not generate enough small primes.\n");
- if (verbose > 1) printf("generated %lu small primes, from 2 to %lu\n", i,
primes[i]);
+ if (verbose > 1) printf("generated %lu small primes, from 2 to %lu\n", i,
(unsigned long)primes[i]);
return primes;
}
@@ -274,14 +276,14 @@ static UV icbrt(UV n)
* this we can avoid caching prime counts and also skip most calls to the
* segment siever.
*/
-static UV bs_prime_count(UV n, UV const* const primes, UV lastidx)
+static UV bs_prime_count(uint32_t n, uint32_t const* const primes, uint32_t
lastidx)
{
UV i, j;
if (n <= 2) return (n == 2);
/* if (n > primes[lastidx]) return _XS_prime_count(2, n); */
if (n >= primes[lastidx]) {
if (n == primes[lastidx]) return lastidx;
- croak("called bspc(%lu) with counts up to %lu\n", n, primes[lastidx]);
+ croak("called bspc(%u) with counts up to %u\n", n, primes[lastidx]);
}
j = lastidx;
if (n < 8480) {
@@ -407,13 +409,13 @@ static void phi_cache_insert(uint32_t x, uint32_t a, IV
sum, cache_t* cache) {
cache->max[a] = cap;
}
if (sum < SHRT_MIN || sum > SHRT_MAX)
- croak("phi(%lu,%lu) 16-bit overflow: sum = %ld\n", x, a, sum);
+ croak("phi(%u,%u) 16-bit overflow: sum = %ld\n", x, a, sum);
if (cache->val[a] == 0)
croak("phi cache allocation failure");
cache->val[a][x] = sum;
}
-static IV _phi3(UV x, UV a, int sign, const UV* const primes, const UV
lastidx, cache_t* cache)
+static IV _phi3(UV x, UV a, int sign, const uint32_t* const primes, const
uint32_t lastidx, cache_t* cache)
{
IV sum;
@@ -612,7 +614,7 @@ static UV phi(UV x, UV a)
UV i, val, sval, lastidx, lastprime;
UV sum = 0;
IV count;
- const UV* primes;
+ const uint32_t* primes;
vcarray_t a1, a2;
vc_t* arr;
cache_t pcache; /* Cache for recursive phi */
@@ -700,10 +702,10 @@ static UV phi(UV x, UV a)
extern UV _XS_meissel_pi(UV n);
/* b = prime_count(isqrt(n)) */
-static UV Pk_2_p(UV n, UV a, UV b, const UV* primes, UV lastprime)
+static UV Pk_2_p(UV n, UV a, UV b, const uint32_t* primes, uint32_t lastidx)
{
UV lastw, lastwpc, i, P2;
- UV lastpc = primes[lastprime];
+ UV lastpc = primes[lastidx];
/* Ensure we have a large enough base sieve */
prime_precalc(isqrt(n / primes[a+1]));
@@ -711,7 +713,7 @@ static UV Pk_2_p(UV n, UV a, UV b, const UV* primes, UV
lastprime)
P2 = lastw = lastwpc = 0;
for (i = b; i > a; i--) {
UV w = n / primes[i];
- lastwpc = (w <= lastpc) ? bs_prime_count(w, primes, lastprime)
+ lastwpc = (w <= lastpc) ? bs_prime_count(w, primes, lastidx)
: lastwpc + _XS_prime_count(lastw+1, w);
lastw = w;
P2 += lastwpc;
@@ -722,7 +724,8 @@ static UV Pk_2_p(UV n, UV a, UV b, const UV* primes, UV
lastprime)
static UV Pk_2(UV n, UV a, UV b)
{
UV lastprime = b*SIEVE_MULT+1;
- const UV* primes = generate_small_primes(lastprime);
+ if (lastprime > 203280221) lastprime = 203280221;
+ const uint32_t* primes = generate_small_primes(lastprime);
UV P2 = Pk_2_p(n, a, b, primes, lastprime);
Safefree(primes);
return P2;
@@ -750,8 +753,8 @@ UV _XS_meissel_pi(UV n)
if (n < SIEVE_LIMIT)
return _XS_prime_count(2, n);
- a = _XS_meissel_pi(icbrt(n)); /* a = floor(n^1/3) */
- b = _XS_meissel_pi(isqrt(n)); /* b = floor(n^1/2) */
+ a = _XS_meissel_pi(icbrt(n)); /* a = Pi(floor(n^1/3)) [max 192725]
*/
+ b = _XS_meissel_pi(isqrt(n)); /* b = Pi(floor(n^1/2)) [max 203280221]
*/
sum = phi(n, a) + a - 1 - Pk_2(n, a, b);
return sum;
@@ -762,7 +765,7 @@ UV _XS_meissel_pi(UV n)
UV _XS_lehmer_pi(UV n)
{
UV z, a, b, c, sum, i, j, lastprime, lastpc, lastw, lastwpc;
- const UV* primes = 0; /* small prime cache, first b=pi(z)=pi(sqrt(n)) */
+ const uint32_t* primes = 0; /* small prime cache, first b=pi(z)=pi(sqrt(n))
*/
DECLARE_TIMING_VARIABLES;
if (n < SIEVE_LIMIT)
@@ -779,9 +782,9 @@ UV _XS_lehmer_pi(UV n)
if (verbose > 0) printf("lehmer %lu stage 1: calculate a,b,c \n", n);
TIMING_START;
z = isqrt(n);
- a = _XS_lehmer_pi(isqrt(z)); /* a = floor(n^1/4) */
- b = _XS_lehmer_pi(z); /* b = floor(n^1/2) */
- c = _XS_lehmer_pi(icbrt(n)); /* c = floor(n^1/3) */
+ a = _XS_lehmer_pi(isqrt(z)); /* a = Pi(floor(n^1/4)) [max 6542]
*/
+ b = _XS_lehmer_pi(z); /* b = Pi(floor(n^1/2)) [max 203280221]
*/
+ c = _XS_lehmer_pi(icbrt(n)); /* c = Pi(floor(n^1/3)) [max 192725]
*/
TIMING_END_PRINT("stage 1")
if (verbose > 0) printf("lehmer %lu stage 2: phi(x,a) (z=%lu a=%lu b=%lu
c=%lu)\n", n, z, a, b, c);
@@ -793,6 +796,7 @@ UV _XS_lehmer_pi(UV n)
* get more than necessary, we can use them to speed up some.
*/
lastprime = b*SIEVE_MULT+1;
+ if (lastprime > 203280221) lastprime = 203280221;
if (verbose > 0) printf("lehmer %lu stage 3: %lu small primes\n", n,
lastprime);
TIMING_START;
primes = generate_small_primes(lastprime);
@@ -840,32 +844,32 @@ UV _XS_lehmer_pi(UV n)
*/
UV _XS_LMO_pi(UV n)
{
- UV n12, n13, a, b, sum, i, j, k, lastprime, P2, S1, S2;
- const UV* primes = 0; /* small prime cache */
+ UV n13, a, b, sum, i, j, k, lastprime, P2, S1, S2;
+ const uint32_t* primes = 0; /* small prime cache */
signed char* mu = 0; /* moebius to n^1/3 */
- UV* lpf = 0; /* least prime factor to n^1/3 */
+ uint32_t* lpf = 0; /* least prime factor to n^1/3 */
cache_t pcache; /* Cache for recursive phi */
DECLARE_TIMING_VARIABLES;
if (n < SIEVE_LIMIT)
return _XS_prime_count(2, n);
- n13 = icbrt(n);
- n12 = isqrt(n);
- a = _XS_lehmer_pi(n13);
- b = _XS_lehmer_pi(n12);
+ n13 = icbrt(n); /* n13 = floor(n^1/3) [max 2642245]
*/
+ a = _XS_lehmer_pi(n13); /* a = Pi(floor(n^1/3)) [max 192725]
*/
+ b = _XS_lehmer_pi(isqrt(n)); /* b = Pi(floor(n^1/2)) [max 203280221]
*/
lastprime = b*SIEVE_MULT+1;
+ if (lastprime > 203280221) lastprime = 203280221;
if (lastprime < n13) lastprime = n13;
primes = generate_small_primes(lastprime);
if (primes == 0) croak("Error generating primes.\n");
New(0, mu, n13+1, signed char);
memset(mu, 1, sizeof(signed char) * (n13+1));
- Newz(0, lpf, n13+1, UV);
+ Newz(0, lpf, n13+1, uint32_t);
mu[0] = 0;
for (i = 1; i <= n13; i++) {
- UV primei = primes[i];
+ uint32_t primei = primes[i];
for (j = primei; j <= n13; j += primei) {
mu[j] = -mu[j];
if (lpf[j] == 0) lpf[j] = primei;
@@ -874,7 +878,7 @@ UV _XS_LMO_pi(UV n)
for (j = k; j <= n13; j += k)
mu[j] = 0;
}
- lpf[1] = UV_MAX; /* Set lpf[1] to max */
+ lpf[1] = 4294967295; /* Set lpf[1] to max */
/* Remove mu[i] == 0 using lpf */
for (i = 1; i <= n13; i++)
@@ -895,7 +899,7 @@ UV _XS_LMO_pi(UV n)
TIMING_START;
for (i = k; i+1 < a; i++) {
- UV p = primes[i+1];
+ uint32_t p = primes[i+1];
/* TODO: #pragma omp parallel for reduction(+: S2) firstprivate(pcache)
schedule(dynamic, 16) */
for (j = (n13/p)+1; j <= n13; j++)
if (lpf[j] > p)
diff --git a/xt/primecount-many.t b/xt/primecount-many.t
new file mode 100644
index 0000000..c93b617
--- /dev/null
+++ b/xt/primecount-many.t
@@ -0,0 +1,145 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/prime_count prime_count_lower prime_count_upper
prime_count_approx/;
+use Digest::SHA qw/sha256_hex/;
+
+my %pivals = (
+ 1000 => 168,
+ 10000 => 1229,
+ 100000 => 9592,
+ 1000000 => 78498,
+ 10000000 => 664579,
+ 100000000 => 5761455,
+ 1000000000 => 50847534,
+ 10000000000 => 455052511,
+ 100000000000 => 4118054813,
+ 1000000000000 => 37607912018,
+ 2000000000000 => 73301896139,
+ 3000000000000 => 108340298703,
+ 4000000000000 => 142966208126,
+ 5000000000000 => 177291661649,
+ 6000000000000 => 211381427039,
+ 7000000000000 => 245277688804,
+ 8000000000000 => 279010070811,
+ 9000000000000 => 312600354108,
+ 10000000000000 => 346065536839,
+ 20000000000000 => 675895909271,
+ 30000000000000 => 1000121668853,
+ 40000000000000 => 1320811971702,
+ 50000000000000 => 1638923764567,
+ 60000000000000 => 1955010428258,
+ 70000000000000 => 2269432871304,
+ 80000000000000 => 2582444113487,
+ 90000000000000 => 2894232250783,
+ 100000000000000 => 3204941750802,
+ 200000000000000 => 6270424651315,
+ 300000000000000 => 9287441600280,
+ 400000000000000 => 12273824155491,
+ 500000000000000 => 15237833654620,
+ 600000000000000 => 18184255291570,
+ 700000000000000 => 21116208911023,
+ 800000000000000 => 24035890368161,
+ 900000000000000 => 26944926466221,
+ 1000000000000000 => 29844570422669,
+ 10000000000000000 => 279238341033925,
+ 20000000000000000 => 547863431950008,
+ 40000000000000000 => 1075292778753150,
+ 100000000000000000 => 2623557157654233,
+ 1000000000000000000 => 24739954287740860,
+ 2000000000000000000 => 48645161281738535,
+ 3000000000000000000 => 72254704797687083,
+ 4000000000000000000 => 95676260903887607,
+ 4185296581467695669 => 100000000000000000,
+ 5000000000000000000 => 118959989688273472,
+ 6000000000000000000 => 142135049412622144,
+ 7000000000000000000 => 165220513980969424,
+ 8000000000000000000 => 188229829247429504,
+ 9000000000000000000 => 211172979243258278,
+10000000000000000000 => 234057667276344607,
+ 524288 => 43390,
+ 1048576 => 82025,
+ 2097152 => 155611,
+ 4194304 => 295947,
+ 8388608 => 564163,
+ 16777216 => 1077871,
+ 33554432 => 2063689,
+ 67108864 => 3957809,
+ 134217728 => 7603553,
+ 268435456 => 14630843,
+ 536870912 => 28192750,
+ 1073741824 => 54400028,
+ 2147483648 => 105097565,
+ 4294967296 => 203280221,
+ 8589934592 => 393615806,
+ 17179869184 => 762939111,
+ 34359738368 => 1480206279,
+ 68719476736 => 2874398515,
+ 137438953472 => 5586502348,
+ 274877906944 => 10866266172,
+ 549755813888 => 21151907950,
+ 1099511627776 => 41203088796,
+ 2199023255552 => 80316571436,
+ 4398046511104 => 156661034233,
+ 8796093022208 => 305761713237,
+ 17592186044416 => 597116381732,
+ 35184372088832 => 1166746786182,
+ 70368744177664 => 2280998753949,
+ 140737488355328 => 4461632979717,
+ 281474976710656 => 8731188863470,
+ 562949953421312 => 17094432576778,
+ 1125899906842624 => 33483379603407,
+ 2251799813685248 => 65612899915304,
+ 4503599627370496 => 128625503610475,
+ 9007199254740992 => 252252704148404,
+ 18014398509481984 => 494890204904784,
+ 36028797018963968 => 971269945245201,
+ 72057594037927936 => 1906879381028850,
+ 144115188075855872 => 3745011184713964,
+ 288230376151711744 => 7357400267843990,
+ 576460752303423488 => 14458792895301660,
+ 1152921504606846976 => 28423094496953330,
+ 2305843009213693952 => 55890484045084135,
+ 4611686018427387904 => 109932807585469973,
+ 9223372036854775808 => 216289611853439384,
+);
+
+plan tests => 5 + 4*scalar(keys %pivals);
+
+# Test prime counts using sampling
+diag "Sampling small prime counts, should take < 1 minute";
+{
+ my $countstr;
+ $countstr = join(" ", map { prime_count($_) } 1 .. 100000);
+ is(sha256_hex($countstr),
"cdbc5c94a927d0d9481cb26b3d3e60c0617a4be65ce9db3075c0363c7a81ef52", "prime
counts 1..10^5");
+ $countstr = join(" ", map { prime_count(100*$_ + ($_%101)) } 1000 .. 100000);
+ is(sha256_hex($countstr),
"73a0b71dedff9611e06fd57e52b88c8afd7f86b5351e4950b2dd5c1d68845b6e", "prime
counts 10^5..10^7 (sample 100)");
+ $countstr = join(" ", map { prime_count(10000*$_ + ($_%9973)) } 1000 ..
10000);
+ is(sha256_hex($countstr),
"d73736c54362136aa0a48bab44b55004b2e63e0d1d03a6cbe1aab42c6a579d0c", "prime
counts 10^7..10^8 (sample 10k)");
+ $countstr = join(" ", map { prime_count(500000*$_ + 250837 + $_) } 200 ..
2000);
+ is(sha256_hex($countstr),
"00a580b2f52b661f065f5ce49bd2aeacb3b169d8903cf824b65731441e40f0b9", "prime
counts 10^8..10^9 (sample 500k)");
+ $countstr = join(" ", map { prime_count(10000000*$_ + 250837 + $_) } 100 ..
1000);
+ is(sha256_hex($countstr),
"9fd78debf4b510ee6d230cabf314ebef5eb253ee63d5df658e45414613f7b8c2", "prime
counts 10^9..10^10 (sample 10M)");
+}
+
+diag "Approximations and limits, should be a few seconds";
+foreach my $n (sort {$a <=> $b} keys %pivals) {
+ my $pin = $pivals{$n};
+ cmp_ok( prime_count_upper($n), '>=', $pin, "Pi($n) <= upper estimate" );
+ cmp_ok( prime_count_lower($n), '<=', $pin, "Pi($n) >= lower estimate" );
+ my $approx = prime_count_approx($n);
+ my $percent_limit = ($n > 1000000000000) ? 0.00005
+ : ($n > 10000000000) ? 0.0002
+ : ($n > 100000000) ? 0.002
+ : ($n > 1000000) ? 0.02
+ : 0.2;
+ cmp_ok( abs($pin - $approx) * (100.0 / $percent_limit), '<=', $pin,
"prime_count_approx($n) within $percent_limit\% of Pi($n)");
+}
+
+diag "Prime counts, will take a very long time";
+foreach my $n (sort {$a <=> $b} keys %pivals) {
+ my $pin = $pivals{$n};
+ is( prime_count($n), $pin, "Pi($n) = $pin" );
+}
--
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
