This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.36 in repository libmath-prime-util-perl.
commit c4ee1370470f26f9752cd12c68c1bc56cb186db1 Author: Dana Jacobsen <d...@acm.org> Date: Sat Jan 11 17:27:39 2014 -0800 Use constants.h. Simplify prev/next prime. Simplify chebyshev functions. Results in smaller object files. --- XS.xs | 6 +-- sieve.c | 17 +++++-- sieve.h | 17 ++++--- util.c | 155 ++++++++++++++++++++-------------------------------------------- util.h | 3 +- 5 files changed, 76 insertions(+), 122 deletions(-) diff --git a/XS.xs b/XS.xs index 7e170d8..6b6a19b 100644 --- a/XS.xs +++ b/XS.xs @@ -777,11 +777,7 @@ _XS_ExponentialIntegral(IN SV* x) } } else { UV uv = SvUV(x); - switch (ix) { - case 4: ret = (NV) _XS_chebyshev_theta(uv); break; - case 5: - default:ret = (NV) _XS_chebyshev_psi(uv); break; - } + ret = (NV) chebyshev_function(uv, ix-4); } RETVAL = ret; OUTPUT: diff --git a/sieve.c b/sieve.c index 8a4afc2..f359747 100644 --- a/sieve.c +++ b/sieve.c @@ -246,7 +246,7 @@ unsigned char* sieve_erat30(UV end) prime = sieve_prefill(mem, 0, max_buf-1); limit = isqrt(end); /* prime*prime can overflow */ - for ( ; prime <= limit; prime = next_prime_in_sieve(mem,prime)) { + for ( ; prime <= limit; prime = next_prime_in_sieve(mem,prime,end)) { UV p2 = prime*prime; UV d = p2 / 30; UV m = p2 - d*30; @@ -293,13 +293,21 @@ unsigned char* sieve_erat30(UV end) int sieve_segment(unsigned char* mem, UV startd, UV endd) { const unsigned char* sieve; - UV limit, slimit, start_base_prime; + UV limit, slimit, start_base_prime, sieve_size; UV startp = 30*startd; UV endp = (endd >= (UV_MAX/30)) ? UV_MAX-2 : 30*endd+29; MPUassert( (mem != 0) && (endd >= startd) && (endp >= startp), "sieve_segment bad arguments"); + /* It's possible we can just use the primary cache */ + sieve_size = get_prime_cache(0, &sieve); + if (sieve_size >= endp) { + memcpy(mem, sieve+startd, endd-startd+1); + release_prime_cache(sieve); + return 1; + } + /* Fill buffer with marked 7, 11, and 13 */ start_base_prime = sieve_prefill(mem, startd, endd); @@ -309,7 +317,10 @@ int sieve_segment(unsigned char* mem, UV startd, UV endd) slimit = limit; if (slimit > BASE_SIEVE_LIMIT) slimit = BASE_SIEVE_LIMIT; /* printf("segment sieve from %"UVuf" to %"UVuf" (aux sieve to %"UVuf")\n", startp, endp, slimit); */ - get_prime_cache(slimit, &sieve); + if (slimit > sieve_size) { + release_prime_cache(sieve); + get_prime_cache(slimit, &sieve); + } START_DO_FOR_EACH_SIEVE_PRIME(sieve, start_base_prime, slimit) { diff --git a/sieve.h b/sieve.h index 70ff0a1..c14cc2d 100644 --- a/sieve.h +++ b/sieve.h @@ -46,18 +46,22 @@ static int is_prime_in_sieve(const unsigned char* sieve, UV p) { #endif #ifdef FUNC_next_prime_in_sieve -/* Warning -- can go off the end of the sieve */ -static UV next_prime_in_sieve(const unsigned char* sieve, UV p) { +/* Will return 0 if it goes past lastp */ +static UV next_prime_in_sieve(const unsigned char* sieve, UV p, UV lastp) { UV d, m; if (p < 7) return (p < 2) ? 2 : (p < 3) ? 3 : (p < 5) ? 5 : 7; d = p/30; m = p - d*30; do { - if (m==29) { d++; m = 1; while (sieve[d] == 0xFF) d++; } - else { m = nextwheel30[m]; } + if (m != 29) { + m = nextwheel30[m]; + } else { + d++; m = 1; + if (d*30 >= lastp) return 0; /* sieves have whole bytes filled */ + } } while (sieve[d] & masktab30[m]); - return(d*30+m); + return d*30+m; } #endif #ifdef FUNC_prev_prime_in_sieve @@ -68,7 +72,8 @@ static UV prev_prime_in_sieve(const unsigned char* sieve, UV p) { d = p/30; m = p - d*30; do { - m = prevwheel30[m]; if (m==29) { if (d == 0) return 0; d--; } + m = prevwheel30[m]; + if (m==29) { if (d == 0) return 0; d--; } } while (sieve[d] & masktab30[m]); return(d*30+m); } diff --git a/util.c b/util.c index 9a3515c..0a51199 100644 --- a/util.c +++ b/util.c @@ -67,6 +67,8 @@ #define FUNC_lcm_ui 1 #define FUNC_ctz 1 #define FUNC_log2floor 1 +#define FUNC_next_prime_in_sieve 1 +#define FUNC_prev_prime_in_sieve 1 #include "util.h" #include "sieve.h" #include "primality.h" @@ -74,6 +76,7 @@ #include "lmo.h" #include "factor.h" #include "mulmod.h" +#include "constants.h" static int _verbose = 0; void _XS_set_verbose(int v) { _verbose = v; } @@ -240,41 +243,20 @@ int _XS_is_prime(UV n) UV _XS_next_prime(UV n) { - UV d, m; + UV d, m, sieve_size, next; const unsigned char* sieve; - UV sieve_size; - - if (n <= 10) { - switch (n) { - case 0: case 1: return 2; break; - case 2: return 3; break; - case 3: case 4: return 5; break; - case 5: case 6: return 7; break; - default: return 11; break; - } - } + if (n < 30*NPRIME_SIEVE30) { - START_DO_FOR_EACH_SIEVE_PRIME(prime_sieve30, n+1, 30*NPRIME_SIEVE30) - return p; - END_DO_FOR_EACH_SIEVE_PRIME; + next = next_prime_in_sieve(prime_sieve30, n,30*NPRIME_SIEVE30); + if (next != 0) return next; } - /* Overflow */ -#if BITS_PER_WORD == 32 - if (n >= UVCONST(4294967291)) return 0; -#else - if (n >= UVCONST(18446744073709551557)) return 0; -#endif + if (n >= MPU_MAX_PRIME) return 0; /* Overflow */ sieve_size = get_prime_cache(0, &sieve); - if (n < sieve_size) { - START_DO_FOR_EACH_SIEVE_PRIME(sieve, n+1, sieve_size) - { release_prime_cache(sieve); return p; } - END_DO_FOR_EACH_SIEVE_PRIME; - /* Not found, so must be larger than the cache size */ - n = sieve_size; - } + next = (n < sieve_size) ? next_prime_in_sieve(sieve, n, sieve_size) : 0; release_prime_cache(sieve); + if (next != 0) return next; d = n/30; m = n - d*30; @@ -292,41 +274,27 @@ UV _XS_next_prime(UV n) UV _XS_prev_prime(UV n) { - UV d, m; const unsigned char* sieve; - UV sieve_size; - - if (n <= 7) - return (n <= 2) ? 0 : (n <= 3) ? 2 : (n <= 5) ? 3 : 5; + UV d, m, prev; - d = n/30; - m = n - d*30; + if (n < 30*NPRIME_SIEVE30) + return prev_prime_in_sieve(prime_sieve30, n); - if (n < 30*NPRIME_SIEVE30) { - /* Don't try to merge this loop with the next one. prime_sieve30 is a - * static, so this can be faster. */ - do { - m = prevwheel30[m]; - if (m==29) d--; - } while (prime_sieve30[d] & masktab30[m]); - return(d*30+m); - } - - sieve_size = get_prime_cache(0, &sieve); - if (n < sieve_size) { - do { - m = prevwheel30[m]; - if (m==29) d--; - } while (sieve[d] & masktab30[m]); + if (n < get_prime_cache(0, &sieve)) { + prev = prev_prime_in_sieve(sieve, n); release_prime_cache(sieve); - } else { - release_prime_cache(sieve); - do { - m = prevwheel30[m]; - if (m==29) d--; - } while (!_is_prime7(d*30+m)); + return prev; } - return(d*30+m); + release_prime_cache(sieve); + + d = n/30; + m = n - d*30; + do { + m = prevwheel30[m]; + if (m==29) d--; + n = d*30+m; + } while (!_is_prime7(n)); + return n; } @@ -678,11 +646,7 @@ static UV _XS_nth_prime_upper(UV n) */ /* 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 + if (n <= MPU_MAX_PRIME_IDX) return MPU_MAX_PRIME; croak("nth_prime_upper(%"UVuf") overflow", n); } @@ -1124,58 +1088,37 @@ UV znlog(UV a, UV g, UV p) { return k; } -long double _XS_chebyshev_theta(UV n) +long double chebyshev_function(UV n, int which) { + long double logp, logn = logl(n); + UV sqrtn = which ? isqrt(n) : 0; /* for theta, p <= sqrtn always false */ KAHAN_INIT(sum); - if (n < 10000) { /* all in one */ - START_DO_FOR_EACH_PRIME(2, n) { - KAHAN_SUM(sum, logl(p)); - } END_DO_FOR_EACH_PRIME - } else { /* Segmented */ - UV seg_base, seg_low, seg_high; - unsigned char* segment; - void* ctx; - KAHAN_SUM(sum, logl(2)); - KAHAN_SUM(sum, logl(3)); - KAHAN_SUM(sum, logl(5)); - ctx = start_segment_primes(7, n, &segment); - while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) { - START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) { - KAHAN_SUM(sum, logl(seg_base + p)); - } END_DO_FOR_EACH_SIEVE_PRIME + if (n < primes_small[NPRIMES_SMALL-1]) { + UV p, pi; + for (pi = 1; (p = primes_small[pi]) <= n; pi++) { + logp = logl(p); + if (p <= sqrtn) logp *= floorl(logn/logp+1e-15); + KAHAN_SUM(sum, logp); } - end_segment_primes(ctx); - } - return sum; -} -long double _XS_chebyshev_psi(UV n) -{ - long double logn = logl(n); - UV sqrtn = isqrt(n); - KAHAN_INIT(sum); - - if (n < 10000) { - START_DO_FOR_EACH_PRIME(2, n) { - long double logp = logl(p); - KAHAN_SUM(sum, (p > (n/p)) ? logp : logp * floorl(logn/logp + 1e-15)); - } END_DO_FOR_EACH_PRIME - return (double) sum; - } - - START_DO_FOR_EACH_PRIME(2, sqrtn) { - long double logp = logl(p); - KAHAN_SUM(sum, logp * floorl(logn/logp + 1e-15)); - } END_DO_FOR_EACH_PRIME - - { /* Segment from sqrt(n) to n */ + } else { UV seg_base, seg_low, seg_high; unsigned char* segment; void* ctx; - ctx = start_segment_primes(sqrtn+1, n, &segment); + if (!which) { + KAHAN_SUM(sum,logl(2)); KAHAN_SUM(sum,logl(3)); KAHAN_SUM(sum,logl(5)); + } else { + KAHAN_SUM(sum, logl(2) * floorl(logn/logl(2) + 1e-15)); + KAHAN_SUM(sum, logl(3) * floorl(logn/logl(3) + 1e-15)); + KAHAN_SUM(sum, logl(5) * floorl(logn/logl(5) + 1e-15)); + } + ctx = start_segment_primes(7, n, &segment); while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) { START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) { - KAHAN_SUM(sum, logl(seg_base + p)); + p += seg_base; + logp = logl(p); + if (p <= sqrtn) logp *= floorl(logn/logp+1e-15); + KAHAN_SUM(sum, logp); } END_DO_FOR_EACH_SIEVE_PRIME } end_segment_primes(ctx); diff --git a/util.h b/util.h index 672f8c7..ed3c8f2 100644 --- a/util.h +++ b/util.h @@ -18,8 +18,7 @@ extern UV _XS_nth_prime(UV x); extern signed char* _moebius_range(UV low, UV high); extern UV* _totient_range(UV low, UV high); extern IV mertens(UV n); -extern long double _XS_chebyshev_theta(UV n); -extern long double _XS_chebyshev_psi(UV n); +extern long double chebyshev_function(UV n, int which); /* 0 = theta, 1 = psi */ extern long double _XS_ExponentialIntegral(long double x); extern long double _XS_LogarithmicIntegral(long double 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