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 <[email protected]> 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 <[email protected]>', - 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 _______________________________________________ Pkg-perl-cvs-commits mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits
