[libmath-prime-util-perl] 165/181: Use constants.h. Simplify prev/next prime. Simplify chebyshev functions. Results in smaller object files.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:18 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list