[libmath-prime-util-perl] 36/50: Fix prime_count bug (added to test suite), add Lehmer prime count
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:38 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.13
in repository libmath-prime-util-perl.
commit c4013aaceca773a5a661aaa600fdb95ea7871f6c
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Nov 15 00:50:09 2012 -0800
Fix prime_count bug (added to test suite), add Lehmer prime count
---
Changes | 6 +-
MANIFEST | 2 +
Makefile.PL | 7 +-
XS.xs | 7 ++
lehmer.c | 257 +++++++++++++++++++++++++++++++++++++++++++++++++
lehmer.h | 9 ++
lib/Math/Prime/Util.pm | 17 ++++
t/13-primecount.t | 2 +
util.c | 16 ++-
9 files changed, 318 insertions(+), 5 deletions(-)
diff --git a/Changes b/Changes
index cf18037..64480bb 100644
--- a/Changes
+++ b/Changes
@@ -1,6 +1,6 @@
Revision history for Perl extension Math::Prime::Util.
-0.12 2 November 2012
+0.12 12 November 2012
- Add bin/primes.pl and bin/factor.pl
@@ -28,6 +28,10 @@ Revision history for Perl extension Math::Prime::Util.
- Complete rewrite of XS p-1 factor routine, includes second stage.
+ - bug fix for prime_count on edge of cache.
+
+ - Add Lehmer prime counting algorithm. This is much faster than sieving.
+
0.11 23 July 2012
- Turn off threading tests on Cygwin, as threads on some Cygwin platforms
give random panics (my Win7 64-bit works fine, XP 32-bit does not).
diff --git a/MANIFEST b/MANIFEST
index eef1b35..608078e 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -14,6 +14,8 @@ cache.h
cache.c
factor.h
factor.c
+lehmer.h
+lehmer.c
sieve.h
sieve.c
util.h
diff --git a/Makefile.PL b/Makefile.PL
index 1592bf7..baa3630 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -7,10 +7,11 @@ WriteMakefile1(
LICENSE => 'perl',
AUTHOR => 'Dana A Jacobsen <dana at acm.org>',
- OBJECT => 'cache.o ' .
- 'factor.o ' .
+ OBJECT => 'cache.o ' .
+ 'factor.o ' .
+ 'lehmer.o ' .
'sieve.o ' .
- 'util.o ' .
+ 'util.o ' .
'XS.o',
LIBS => ['-lm'],
diff --git a/XS.xs b/XS.xs
index 25dadce..84b0066 100644
--- a/XS.xs
+++ b/XS.xs
@@ -8,6 +8,7 @@
#include "sieve.h"
#include "util.h"
#include "factor.h"
+#include "lehmer.h"
MODULE = Math::Prime::Util PACKAGE = Math::Prime::Util
@@ -40,6 +41,12 @@ _XS_prime_count(IN UV low, IN UV high = 0)
RETVAL
UV
+_XS_legendre_pi(IN UV n)
+
+UV
+_XS_lehmer_pi(IN UV n)
+
+UV
_XS_nth_prime(IN UV n)
int
diff --git a/lehmer.c b/lehmer.c
new file mode 100644
index 0000000..23844bd
--- /dev/null
+++ b/lehmer.c
@@ -0,0 +1,257 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "lehmer.h"
+#include "util.h"
+#include "cache.h"
+#include "sieve.h"
+
+static int const verbose = 0;
+
+static IV mapes(IV x, UV a)
+{
+ IV val;
+ if (a == 0) return x;
+ if (a == 1) return x-x/2;
+ val = x-x/2-x/3+x/6;
+ if (a >= 3) val += -x/5+x/10+x/15-x/30;
+ if (a >= 4) val += -x/7+x/14+x/21-x/42+x/35-x/70-x/105+x/210;
+ if (a >= 5) val += -x/11+x/22+x/33-x/66+x/55-x/110-x/165+x/330+x/77-x/154-x/231+x/462-x/385+x/770+x/1155-x/2310;
+ if (a >= 6) val += -x/13+x/26+x/39-x/78+x/65-x/130-x/195+x/390+x/91-x/182-x/273+x/546-x/455+x/910+x/1365-x/2730+x/143-x/286-x/429+x/858-x/715+x/1430+x/2145-x/4290-x/1001+x/2002+x/3003-x/6006+x/5005-x/10010-x/15015+x/30030;
+ if (a >= 7) val += -x/17+x/34+x/51-x/102+x/85-x/170-x/255+x/510+x/119-x/238-x/357+x/714-x/595+x/1190+x/1785-x/3570+x/187-x/374-x/561+x/1122-x/935+x/1870+x/2805-x/5610-x/1309+x/2618+x/3927-x/7854+x/6545-x/13090-x/19635+x/39270+x/221-x/442-x/663+x/1326-x/1105+x/2210+x/3315-x/6630-x/1547+x/3094+x/4641-x/9282+x/7735-x/15470-x/23205+x/46410-x/2431+x/4862+x/7293-x/14586+x/12155-x/24310-x/36465+x/72930+x/17017-x/34034-x/51051+x/102102-x/85085+x/170170+x/255255-x/510510;
+ return val;
+}
+
+/* heap of values and counts, stored by value */
+typedef struct {
+ UV v;
+ IV c;
+} vc_t;
+
+typedef struct {
+ vc_t* array;
+ UV array_size;
+ UV N;
+} heap_t;
+
+static void heap_insert(heap_t* h, vc_t v)
+{
+ UV k;
+ UV n = ++(h->N);
+ vc_t* a = h->array;
+ if (n >= h->array_size) {
+ UV new_size = (h->array_size == 0) ? 20000 : 2 * h->array_size;
+ if (verbose>2) printf("REALLOCing %p, new size %lu\n", h->array, new_size);
+ h->array = realloc( h->array, new_size * sizeof(vc_t) );
+ if (h->array == 0)
+ croak("could not allocate heap\n");
+ a = h->array;
+ h->array_size = new_size-1;
+ a[0].v = UV_MAX; a[0].c = 0;
+ }
+ a[n] = v;
+ k = n;
+ while (a[k/2].v <= v.v) { /* upheap */
+ a[k] = a[k/2];
+ k /= 2;
+ }
+ a[k] = v;
+}
+
+static vc_t heap_remove(heap_t* h)
+{
+ UV k;
+ UV n = h->N;
+ vc_t* a = h->array;
+ vc_t v, top;
+ if (n == 0)
+ croak("removing from empty heap\n");
+ top = a[1];
+ v = a[1] = a[n];
+ n = --(h->N); /* downheap */
+ k = 1;
+ while (k <= n/2) {
+ UV j = k+k;
+ if (j < n && a[j].v < a[j+1].v) j++;
+ if (v.v >= a[j].v) break;
+ a[k] = a[j];
+ k = j;
+ }
+ a[k] = v;
+ return top;
+}
+static vc_t heap_remove_coalesce(heap_t* h)
+{
+ vc_t v = heap_remove(h);
+ /* get rest of entries with same value */
+ while (h->N > 0 && h->array[1].v == v.v) {
+ v.c += h->array[1].c;
+ (void) heap_remove(h);
+ }
+ return v;
+}
+static heap_t heap_create(void)
+{
+ heap_t h;
+ h.array = 0;
+ h.array_size = 0;
+ h.N = 0;
+ return h;
+}
+static void heap_destroy(heap_t* h)
+{
+ if (h->array != 0) {
+ if (verbose > 2) printf("FREE %p\n", h->array);
+ free(h->array);
+ }
+ h->array = 0;
+ h->array_size = 0;
+ h->N = 0;
+}
+static void heap_empty(heap_t* h)
+{
+ h->N = 0;
+}
+
+
+static UV phi(UV x, UV a)
+{
+ UV i;
+ IV primea;
+ IV sum = 0;
+ if (x < _XS_nth_prime(a+1)) return (x > 0) ? 1 : 0;
+ if (a == 1) return ((x+1)/2);
+ if (a <= 7) return mapes(x, a);
+
+ primea = _XS_nth_prime(a);
+
+ heap_t h1 = heap_create();
+ heap_t h2 = heap_create();
+ vc_t v;
+
+ v.v = x; v.c = 1;
+ heap_insert(&h1, v);
+ while (a > 5) {
+ heap_empty(&h2);
+ /* Walk heap */
+ while (h1.N > 0) {
+ UV sval;
+ v = heap_remove_coalesce(&h1);
+ if (v.c == 0)
+ continue;
+ heap_insert(&h2, v);
+ sval = v.v / primea;
+ if (sval >= primea) {
+ v.v = sval;
+ v.c = -v.c;
+ heap_insert(&h2, v);
+ } else {
+ sum -= v.c;
+ }
+ }
+ { heap_t t = h1; h1 = h2; h2 = t; }
+ a--;
+ primea = _XS_prev_prime(primea);
+ }
+ heap_destroy(&h2);
+ while (h1.N > 0) {
+ v = heap_remove_coalesce(&h1);
+ if (v.c != 0)
+ sum += v.c * mapes(v.v, a);
+ }
+ heap_destroy(&h1);
+ return sum;
+}
+
+UV _XS_legendre_pi(UV n)
+{
+ if (n < 30000)
+ return _XS_prime_count(2, n);
+
+ /* Legendre */
+ UV a = _XS_legendre_pi(sqrt(n));
+ prime_precalc(a);
+ return phi(n, a) + a - 1;
+}
+
+UV _XS_lehmer_pi(UV n)
+{
+ UV z, a, b, c, sum, i, j, lastw, lastwpc;
+ UV* primea = 0; /* small prime cache, first b=pi(z)=pi(sqrt(n)) primes */
+ UV* pca = 0; /* small prime count cache, first z=sqrt(n) numbers */
+ if (n < 1000000)
+ return _XS_prime_count(2, n);
+
+ if (verbose > 0) printf("lehmer stage 1 - %lu\n", n);
+ z = sqrt((double)n+0.5);
+ a = _XS_lehmer_pi(sqrt((double)z)+0.5);
+ b = _XS_lehmer_pi(z);
+ c = _XS_lehmer_pi(pow((double)n, 1.0/3.0)+0.5);
+ if (verbose > 0) printf("lehmer stage 2 - %lu (z=%lu a=%lu b=%lu c=%lu)\n", n, z, a, b, c);
+ sum = phi(n, a) + ( (b+a-2) * (b-a+1) / 2);
+
+ if (verbose > 0) printf("lehmer stage 3 - n=%lu a=%lu\n", n, a);
+ /* Calculate the first b primes and first z primecounts */
+ {
+ UV count = 0;
+ const unsigned char* sieve;
+ if (get_prime_cache(z, &sieve) < z) {
+ release_prime_cache(sieve);
+ croak("Could not generate sieve for %"UVuf, z);
+ }
+ primea = (UV*) malloc( (b+4) * sizeof(UV) );
+ if (primea == 0)
+ croak("Can not allocate small primes\n");
+ primea[0] = 0;
+ primea[1] = 2;
+ primea[2] = 3;
+ primea[3] = 5;
+ i = 3;
+ START_DO_FOR_EACH_SIEVE_PRIME( sieve, 7, z ) {
+ if (i > b) break;
+ primea[++i] = p;
+ } END_DO_FOR_EACH_SIEVE_PRIME
+ release_prime_cache(sieve);
+ if (i < b)
+ croak("Did not generate enough small primes. Bug in Lehmer code.\n");
+ if (verbose > 1) printf("generated first %lu small primes, from 2 to %lu\n", i, z);
+
+ pca = (UV*) malloc( (z+4) * sizeof(UV) );
+ if (pca == 0)
+ croak("Can not allocate small prime counts\n");
+ for (i = 0; i <= z && count+1 <= b; i++) {
+ if (i >= primea[count+1])
+ count++;
+ pca[i] = count;
+ }
+ if (verbose > 1) printf("generated first %lu prime counts, from 2 to %lu\n", count, z);
+ }
+
+ /* Ensure we have the base sieve for prime_count ( n/primea[i] ) */
+ prime_precalc( sqrt( n / primea[a+1] ) );
+
+
+ if (verbose > 0) printf("lehmer stage 4 - n=%lu a=%lu\n", n, a);
+ /* Reverse the i loop so w increases. Count w in segments. */
+ lastw = 0;
+ lastwpc = 0;
+ for (i = b; i >= a+1; i--) {
+ UV w = n / primea[i];
+ lastwpc += _XS_prime_count(lastw+1, w);
+ lastw = w;
+ sum = sum - lastwpc;
+ if (i <= c) {
+ UV bi = pca[(UV) (sqrt(w) + 0.5)];
+ for (j = i; j <= bi; j++) {
+ sum = sum - pca[w / primea[j]] + j - 1;
+ }
+ }
+ }
+ if (pca != 0)
+ free(pca);
+ if (primea != 0)
+ free(primea);
+ return sum;
+}
diff --git a/lehmer.h b/lehmer.h
new file mode 100644
index 0000000..410feb0
--- /dev/null
+++ b/lehmer.h
@@ -0,0 +1,9 @@
+#ifndef MPU_LEHMER_H
+#define MPU_LEHMER_H
+
+#include "ptypes.h"
+
+extern UV _XS_legendre_pi(UV n);
+extern UV _XS_lehmer_pi(UV n);
+
+#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 606b487..2e0e7e0 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -20,6 +20,7 @@ our @EXPORT_OK = qw(
primes
next_prime prev_prime
prime_count prime_count_lower prime_count_upper prime_count_approx
+ prime_count_legendre prime_count_lehmer
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime random_nbit_prime random_maurer_prime
primorial pn_primorial
@@ -867,6 +868,22 @@ sub prime_count {
return Math::Prime::Util::PP::prime_count($low,$high);
}
+sub prime_count_legendre {
+ my($n) = @_;
+ return 0 if $n <= 0;
+ _validate_positive_integer($n);
+
+ return _XS_legendre_pi($n) if $n <= $_XS_MAXVAL;
+}
+
+sub prime_count_lehmer {
+ my($n) = @_;
+ return 0 if $n <= 0;
+ _validate_positive_integer($n);
+
+ return _XS_lehmer_pi($n) if $n <= $_XS_MAXVAL;
+}
+
sub nth_prime {
my($n) = @_;
_validate_positive_integer($n);
diff --git a/t/13-primecount.t b/t/13-primecount.t
index 475cea7..23765db 100644
--- a/t/13-primecount.t
+++ b/t/13-primecount.t
@@ -21,6 +21,8 @@ my %pivals32 = (
10000000 => 664579,
100000000 => 5761455,
1000000000 => 50847534,
+ 30239 => 3269,
+ 30249 => 3270,
60067 => 6062,
65535 => 6542,
16777215 => 1077871,
diff --git a/util.c b/util.c
index 766274d..96e07cc 100644
--- a/util.c
+++ b/util.c
@@ -427,6 +427,20 @@ UV _XS_prime_count(UV low, UV high)
if (low > high) return count;
+ if (low == 7) {
+ if (high > (1UL << 42)) { count += 156661034233-3; low = 1UL<<42; }
+ else if (high > (1UL << 39)) { count += 21151907950-3; low = 1UL<<39; }
+ else if (high > (1UL << 36)) { count += 2874398515-3; low = 1UL<<36; }
+ else if (high > (1UL << 33)) { count += 393615806-3; low = 1UL<<33; }
+ else if (high > (1UL << 30)) { count += 54400028-3; low = 1UL<<30; }
+ else if (high > (1UL << 27)) { count += 7603553-3; low = 1UL<<27; }
+ else if (high > (1UL << 24)) { count += 1077871-3; low = 1UL<<24; }
+ else if (high > (1UL << 20)) { count += 82025-3; low = 1UL<<20; }
+ else if (high > (1UL << 18)) { count += 23000-3; low = 1UL<<18; }
+ else if (high > (1UL << 16)) { count += 6542-3; low = 1UL<<16; }
+ else if (high > (1UL << 14)) { count += 1900-3; low = 1UL<<14; }
+ }
+
low_d = low/30;
high_d = high/30;
@@ -443,7 +457,7 @@ UV _XS_prime_count(UV low, UV high)
/* Count all the primes in the primary cache in our range */
count += count_segment_ranged(cache_sieve, segment_size, low, high);
- if (high_d <= segment_size) {
+ if (high_d < segment_size) {
release_prime_cache(cache_sieve);
return count;
}
--
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