This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.25
in repository libmath-prime-util-perl.
commit 718c44d785f1557512e786e4c235957a9f9203cb
Author: Dana Jacobsen <d...@acm.org>
Date: Mon Mar 11 18:31:36 2013 -0700
New internal macro to loop a..b using primary sieve
---
Changes | 3 +++
TODO | 7 -------
XS.xs | 16 +++-------------
factor.c | 41 +++++++++++++++++++----------------------
lehmer.c | 15 ++++-----------
sieve.h | 35 +++++++++++++++++++++++++++++++++++
util.c | 61 ++++++++++++++-----------------------------------------------
7 files changed, 78 insertions(+), 100 deletions(-)
diff --git a/Changes b/Changes
index 71c04c9..fdc9b74 100644
--- a/Changes
+++ b/Changes
@@ -5,6 +5,9 @@ Revision history for Perl extension Math::Prime::Util.
- Speed up p-1 stage 2 factoring. This makes factor() about 20% faster
for 19 digit semiprimes.
+ - New internal macro to loop over primary sieve starting at 2. Simplifies
+ code in quite a few places.
+
0.24 10 March 2013
- Fix compilation with old pre-C99 strict compilers (decl after statement).
diff --git a/TODO b/TODO
index dc0ddd7..a75cc17 100644
--- a/TODO
+++ b/TODO
@@ -41,11 +41,4 @@
- More efficient totient segment. Do we really need primes to n/2?
-- Add a more complex system for allowing something like:
- FOREACH_PRIME(low, high) {
- ....
- } END_FOREACH_PRIME
- that (1) handles any low / high, and (2) segments. Take segment_primes
- from XS.xs to do this.
-
- Implement S2 calculation for LMO prime count.
diff --git a/XS.xs b/XS.xs
index 5f2f42d..ff46ecc 100644
--- a/XS.xs
+++ b/XS.xs
@@ -102,19 +102,9 @@ sieve_primes(IN UV low, IN UV high)
AV* av = newAV();
CODE:
if (low <= high) {
- if (get_prime_cache(high, &sieve) < high) {
- release_prime_cache(sieve);
- croak("Could not generate sieve for %"UVuf, high);
- } else {
- if ((low <= 2) && (high >= 2)) { av_push(av, newSVuv( 2 )); }
- if ((low <= 3) && (high >= 3)) { av_push(av, newSVuv( 3 )); }
- if ((low <= 5) && (high >= 5)) { av_push(av, newSVuv( 5 )); }
- if (low < 7) { low = 7; }
- START_DO_FOR_EACH_SIEVE_PRIME( sieve, low, high ) {
- av_push(av,newSVuv(p));
- } END_DO_FOR_EACH_SIEVE_PRIME
- release_prime_cache(sieve);
- }
+ START_DO_FOR_EACH_PRIME(low, high) {
+ av_push(av,newSVuv(p));
+ } END_DO_FOR_EACH_PRIME
}
RETVAL = newRV_noinc( (SV*) av );
OUTPUT:
diff --git a/factor.c b/factor.c
index 75e496b..26622cf 100644
--- a/factor.c
+++ b/factor.c
@@ -66,7 +66,7 @@ int factor(UV n, UV *factors)
}
/* This p-1 gets about 2/3 of what makes it through the above */
if (!split_success) {
- split_success = pminus1_factor(n, tofac_stack+ntofac, 4000, 80000)-1;
+ split_success = pminus1_factor(n, tofac_stack+ntofac, 5000, 80000)-1;
if (verbose) printf("pminus1 %d\n", split_success);
}
/* Some rounds of HOLF, good for close to perfect squares */
@@ -632,19 +632,21 @@ int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
UV sqrtB1 = isqrt(B1);
MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in pminus1_factor");
- for (q = 2; q <= sqrtB1; q = _XS_next_prime(q)) {
- UV k = q*q;
- UV kmin = B1/q;
+ START_DO_FOR_EACH_PRIME(2, sqrtB1) {
+ UV k = p*p;
+ UV kmin = B1/p;
while (k <= kmin)
- k *= q;
+ k *= p;
a = powmod(a, k, n);
- }
+ q = p;
+ } END_DO_FOR_EACH_PRIME
if (a == 0) { factors[0] = n; return 1; }
f = gcd_ui(a-1, n);
if (f == 1) {
savea = a;
saveq = q;
- for (; q <= B1; q = _XS_next_prime(q)) {
+ START_DO_FOR_EACH_PRIME(q+1, B1) {
+ q = p;
a = powmod(a, q, n);
if ( (j++ % 32) == 0) {
if (a == 0 || gcd_ui(a-1, n) != 1)
@@ -652,23 +654,24 @@ int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
savea = a;
saveq = q;
}
- }
+ } END_DO_FOR_EACH_PRIME
if (a == 0) { factors[0] = n; return 1; }
f = gcd_ui(a-1, n);
}
/* If we found more than one factor in stage 1, backup and single step */
if (f == n) {
a = savea;
- for (q = saveq; q <= B1; q = _XS_next_prime(q)) {
- UV k = q;
- UV kmin = B1/q;
+ START_DO_FOR_EACH_PRIME(saveq, B1) {
+ UV k = p;
+ UV kmin = B1/p;
while (k <= kmin)
- k *= q;
+ k *= p;
a = powmod(a, k, n);
f = gcd_ui(a-1, n);
+ q = p;
if (f != 1)
break;
- }
+ } END_DO_FOR_EACH_PRIME
/* If f == n again, we could do:
* for (savea = 3; f == n && savea < 100; savea = _XS_next_prime(savea)) {
* a = savea;
@@ -681,8 +684,7 @@ int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
}
/* STAGE 2 */
- if (f == 1 && B2 > B1 && q >= 5) {
- const unsigned char* sieve;
+ if (f == 1 && B2 > B1) {
UV bm = a;
UV b = 1;
UV bmdiff;
@@ -699,11 +701,7 @@ int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
a = powmod(a, q, n);
j = 1;
- if (get_prime_cache(B2, &sieve) < B2) {
- release_prime_cache(sieve);
- croak("Could not generate sieve for %"UVuf, B2);
- }
- START_DO_FOR_EACH_SIEVE_PRIME( sieve, q+1, B2 ) {
+ START_DO_FOR_EACH_PRIME( q+1, B2 ) {
UV lastq = q;
UV qdiff;
q = p;
@@ -729,8 +727,7 @@ int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
if (f != 1)
break;
}
- } END_DO_FOR_EACH_SIEVE_PRIME
- release_prime_cache(sieve);
+ } END_DO_FOR_EACH_PRIME
f = gcd_ui(b, n);
}
if ( (f != 1) && (f != n) ) {
diff --git a/lehmer.c b/lehmer.c
index 267787b..1417e2e 100644
--- a/lehmer.c
+++ b/lehmer.c
@@ -183,9 +183,8 @@ static UV* generate_small_primes(UV n)
* Remember to free when done. */
static UV* generate_small_primes(UV n)
{
- const unsigned char* sieve;
UV* primes;
- UV i;
+ UV i = 0;
double fn = (double)n;
double flogn = log(fn);
double flog2n = log(flogn);
@@ -195,20 +194,14 @@ static UV* generate_small_primes(UV n)
(n >= 18) ? (UV) ceil(fn*(flogn+flog2n-1.0+((flog2n+0.30)/flogn)))
: 59;
- if (get_prime_cache(nth_prime, &sieve) < nth_prime) {
- release_prime_cache(sieve);
- croak("Could not generate sieve for %"UVuf, nth_prime);
- }
New(0, primes, n+1, UV);
if (primes == 0)
croak("Can not allocate small primes\n");
- primes[0] = 0; primes[1] = 2; primes[2] = 3; primes[3] = 5;
- i = 3;
- START_DO_FOR_EACH_SIEVE_PRIME( sieve, 7, nth_prime ) {
+ primes[0] = 0;
+ START_DO_FOR_EACH_PRIME(2, nth_prime) {
if (i >= n) break;
primes[++i] = p;
- } END_DO_FOR_EACH_SIEVE_PRIME
- release_prime_cache(sieve);
+ } 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]);
diff --git a/sieve.h b/sieve.h
index 8e1f3f9..6dca024 100644
--- a/sieve.h
+++ b/sieve.h
@@ -81,5 +81,40 @@ static UV prev_prime_in_sieve(const unsigned char* sieve, UV
p) {
} \
}
+#define START_DO_FOR_EACH_PRIME(a, b) \
+ { \
+ const unsigned char* sieve_; \
+ UV d_ = 0; \
+ UV m_ = 7; \
+ UV p = a; \
+ UV l_ = b; \
+ if (get_prime_cache(l_, &sieve_) < l_) { \
+ release_prime_cache(sieve_); \
+ croak("Could not generate sieve for %"UVuf, l_); \
+ } \
+ if (p <= 5) { \
+ p = (p <= 2) ? 2 : (p <= 3) ? 3 : 5; \
+ } else { \
+ d_ = p/30; \
+ m_ = p-d_*30; \
+ m_ += distancewheel30[m_]; \
+ p = d_*30 + m_; \
+ } \
+ while ( p <= l_ ) { \
+ if (p < 7 || !(sieve_[d_] & masktab30[m_]))
+
+#define RETURN_FROM_EACH_PRIME(x) \
+ do { release_prime_cache(sieve_); return(x); } while (0)
+
+#define END_DO_FOR_EACH_PRIME \
+ if (p < 7) { \
+ p += 1 + (p > 2); \
+ } else { \
+ m_ = nextwheel30[m_]; if (m_ == 1) { d_++; } \
+ p = d_*30+m_; \
+ } \
+ } \
+ release_prime_cache(sieve_); \
+ }
#endif
diff --git a/util.c b/util.c
index 3eb4929..c7c3505 100644
--- a/util.c
+++ b/util.c
@@ -769,7 +769,7 @@ UV _XS_nth_prime(UV n)
char* _moebius_range(UV lo, UV hi)
{
char* mu;
- UV i, p;
+ UV i;
UV sqrtn = isqrt(hi);
#if 1
IV* A;
@@ -786,14 +786,13 @@ char* _moebius_range(UV lo, UV hi)
croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1);
for (i = lo; i <= hi; i++)
A[i-lo] = 1;
- prime_precalc(sqrtn);
- for (p = 2; p <= sqrtn; p = _XS_next_prime(p)) {
+ START_DO_FOR_EACH_PRIME(2, sqrtn) {
UV p2 = p*p;
for (i = PGTLO(p2, lo); i <= hi; i += p2)
A[i-lo] = 0;
for (i = PGTLO(p, lo); i <= hi; i += p)
A[i-lo] *= -(IV)p;
- }
+ } END_DO_FOR_EACH_PRIME
New(0, mu, hi-lo+1, char);
if (mu == 0)
croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1);
@@ -807,6 +806,7 @@ char* _moebius_range(UV lo, UV hi)
Safefree(A);
#endif
#if 0
+ UV p;
/* Simple char method, Needs way too many primes */
New(0, mu, hi-lo+1, char);
if (mu == 0)
@@ -827,6 +827,7 @@ char* _moebius_range(UV lo, UV hi)
* (1) I'm using the log function, which should be fixed (easy).
* (2) it doesn't work. try 64-101 vs. 64 vs. 100.
*/
+ UV p;
unsigned char* A;
New(0, A, hi-lo+1, unsigned char);
if (A == 0)
@@ -862,7 +863,6 @@ char* _moebius_range(UV lo, UV hi)
UV* _totient_range(UV lo, UV hi) {
UV* totients;
- const unsigned char* sieve;
UV i, sievehi;
if (hi < lo) croak("_totient_range error hi %lu < lo %lu\n", hi, lo);
New(0, totients, hi-lo+1, UV);
@@ -872,18 +872,10 @@ UV* _totient_range(UV lo, UV hi) {
totients[i-lo] = i;
sievehi = hi/2;
for (i=P2GTLO(2*2,2,lo); i <= hi; i += 2) totients[i-lo] -= totients[i-lo]/2;
- for (i=P2GTLO(2*3,3,lo); i <= hi; i += 3) totients[i-lo] -= totients[i-lo]/3;
- for (i=P2GTLO(2*5,5,lo); i <= hi; i += 5) totients[i-lo] -= totients[i-lo]/5;
- if (get_prime_cache(sievehi, &sieve) < sievehi) {
- release_prime_cache(sieve);
- croak("Could not generate sieve for %"UVuf, sievehi);
- } else {
- START_DO_FOR_EACH_SIEVE_PRIME( sieve, 7, sievehi ) {
- for (i = P2GTLO(2*p,p,lo); i <= hi; i += p)
- totients[i-lo] -= totients[i-lo]/p;
- } END_DO_FOR_EACH_SIEVE_PRIME
- release_prime_cache(sieve);
- }
+ START_DO_FOR_EACH_PRIME(3, sievehi) {
+ for (i = P2GTLO(2*p,p,lo); i <= hi; i += p)
+ totients[i-lo] -= totients[i-lo]/p;
+ } END_DO_FOR_EACH_PRIME
for (i = lo; i <= hi; i++)
if (totients[i-lo] == i)
totients[i-lo] = i-1;
@@ -951,49 +943,24 @@ IV _XS_mertens(UV n) {
double _XS_chebyshev_theta(UV n)
{
- const unsigned char* sieve;
KAHAN_INIT(sum);
-
- if (n >= 2) KAHAN_SUM(sum, 0.6931471805599453094172321214581765680755L);
- if (n >= 3) KAHAN_SUM(sum, 1.0986122886681096913952452369225257046475L);
- if (n >= 5) KAHAN_SUM(sum, 1.6094379124341003746007593332261876395256L);
- if (n < 7) return (double) sum;
-
- if (get_prime_cache(n, &sieve) < n) {
- release_prime_cache(sieve);
- croak("Could not generate sieve for %"UVuf, n);
- }
- START_DO_FOR_EACH_SIEVE_PRIME(sieve, 7, n) {
+ START_DO_FOR_EACH_PRIME(2, n) {
KAHAN_SUM(sum, logl(p));
- } END_DO_FOR_EACH_SIEVE_PRIME
- release_prime_cache(sieve);
+ } END_DO_FOR_EACH_PRIME
return (double) sum;
}
double _XS_chebyshev_psi(UV n)
{
- const unsigned char* sieve;
- UV prime, mults_are_one;
+ UV mults_are_one = 0;
long double logn, logp;
KAHAN_INIT(sum);
logn = logl(n);
- for (prime = 2; prime <= 5 && prime <= n; prime += prime-1) {
- logp = logl(prime);
- KAHAN_SUM(sum, logp * floorl(logn/logp + 1e-15));
- }
- if (n < 7) return (double) sum;
-
- if (get_prime_cache(n, &sieve) < n) {
- release_prime_cache(sieve);
- croak("Could not generate sieve for %"UVuf, n);
- }
- mults_are_one = 0;
- START_DO_FOR_EACH_SIEVE_PRIME(sieve, 7, n) {
+ START_DO_FOR_EACH_PRIME(2, n) {
logp = logl(p);
if (!mults_are_one && p > (n/p)) mults_are_one = 1;
KAHAN_SUM(sum, (mults_are_one) ? logp : logp * floorl(logn/logp + 1e-15));
- } END_DO_FOR_EACH_SIEVE_PRIME
- release_prime_cache(sieve);
+ } END_DO_FOR_EACH_PRIME
return (double) sum;
}
--
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
