[libmath-prime-util-perl] 169/181: remove _is_prime7, remove _XS_ header from a couple C functions
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 64fe8219f1127a7de42241838a5ba8c14c132b81
Author: Dana Jacobsen <dana at acm.org>
Date: Sun Jan 12 05:46:27 2014 -0800
remove _is_prime7, remove _XS_ header from a couple C functions
---
XS.xs | 16 ++++++-------
aks.c | 2 +-
factor.c | 4 ++--
lmo.c | 4 ++--
mulmod.h | 5 ++--
primality.c | 76 +++++++++++++++++++++----------------------------------------
primality.h | 2 +-
sieve.c | 4 ++--
util.c | 62 +++++++++++++++----------------------------------
util.h | 12 ++--------
10 files changed, 65 insertions(+), 122 deletions(-)
diff --git a/XS.xs b/XS.xs
index 6b6a19b..d1ba176 100644
--- a/XS.xs
+++ b/XS.xs
@@ -379,9 +379,9 @@ sieve_primes(IN UV low, IN UV high)
av_push(av,newSVuv(p));
} END_DO_FOR_EACH_PRIME
} else if (ix == 1) { /* Trial */
- for (low = _XS_next_prime(low-1);
+ for (low = next_prime(low-1);
low <= high && low != 0;
- low = _XS_next_prime(low) ) {
+ low = next_prime(low) ) {
av_push(av,newSVuv(low));
}
} else if (ix == 2) { /* Erat with private memory */
@@ -596,8 +596,8 @@ next_prime(IN SV* svn)
} else {
UV ret;
switch (ix) {
- case 0: ret = _XS_next_prime(n); break;
- case 1: ret = (n < 3) ? 0 : _XS_prev_prime(n); break;
+ case 0: ret = next_prime(n); break;
+ case 1: ret = (n < 3) ? 0 : prev_prime(n); break;
case 2:
default:ret = _XS_nth_prime(n); break;
}
@@ -960,7 +960,7 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
(beg >= UVCONST( 1000000000000) && end-beg < 17000) ||
#endif
((end-beg) < 500) ) { /* MULTICALL next prime */
- for (beg = _XS_next_prime(beg-1); beg <= end && beg != 0; beg = _XS_next_prime(beg)) {
+ for (beg = next_prime(beg-1); beg <= end && beg != 0; beg = next_prime(beg)) {
sv_setuv(svarg, beg);
MULTICALL;
}
@@ -1036,8 +1036,8 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
}
/* Find the two primes that bound their interval. */
/* If beg or end are >= max_prime, then this will die. */
- prevprime = _XS_prev_prime(beg);
- nextprime = _XS_next_prime(end);
+ prevprime = prev_prime(beg);
+ nextprime = next_prime(end);
ctx = start_segment_primes(beg, nextprime, &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 ) {
@@ -1059,7 +1059,7 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
if (beg <= end) {
beg = (beg <= 4) ? 3 : beg-1;
while (beg++ < end) {
- if (!_XS_is_prob_prime(beg)) {
+ if (!is_prob_prime(beg)) {
sv_setuv(svarg, beg);
PUSHMARK(SP);
call_sv((SV*)cv, G_VOID|G_DISCARD);
diff --git a/aks.c b/aks.c
index 5ab49f5..00e8d29 100644
--- a/aks.c
+++ b/aks.c
@@ -50,7 +50,7 @@ static int is_perfect_power(UV n) {
b = isqrt(n);
if (b*b == n) return 1; /* perfect square */
last = log2floor(n-1) + 1;
- for (b = 3; b < last; b = _XS_next_prime(b)) {
+ for (b = 3; b < last; b = next_prime(b)) {
UV root = (UV) (pow(n, 1.0 / (double)b) + 0.5);
if ( ((UV)(pow(root, b)+0.5)) == n) return 1;
}
diff --git a/factor.c b/factor.c
index 8208385..d1259b1 100644
--- a/factor.c
+++ b/factor.c
@@ -593,9 +593,9 @@ int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
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)) {
+ * for (savea = 3; f == n && savea < 100; savea = next_prime(savea)) {
* a = savea;
- * for (q = 2; q <= B1; q = _XS_next_prime(q)) {
+ * for (q = 2; q <= B1; q = next_prime(q)) {
* ...
* }
* }
diff --git a/lmo.c b/lmo.c
index eed5cf4..3d40af5 100644
--- a/lmo.c
+++ b/lmo.c
@@ -297,12 +297,12 @@ static UV _phi_recurse(UV x, UV a) {
UV pa = _XS_nth_prime(a);
for (i = c+1; i <= a; i++) {
UV xp;
- p = _XS_next_prime(p);
+ p = next_prime(p);
xp = x/p;
if (xp < p) {
while (x < pa) {
a--;
- pa = _XS_prev_prime(pa);
+ pa = prev_prime(pa);
}
return (sum - a + i - 1);
}
diff --git a/mulmod.h b/mulmod.h
index da56608..b826e98 100644
--- a/mulmod.h
+++ b/mulmod.h
@@ -82,6 +82,7 @@
UV r = 0;
if (a >= n) a %= n; /* Careful attention from the caller should make */
if (b >= n) b %= n; /* these unnecessary. */
+ if ((a|b) < HALF_WORD) return (a*b) % n;
if (a < b) { UV t = a; a = b; b = t; }
if (n <= (UV_MAX>>1)) {
while (b > 0) {
@@ -99,8 +100,8 @@
return r;
}
- #define mulmod(a,b,n) ((((a)|(b)) < HALF_WORD) ? ((a)*(b))%(n):_mulmod(a,b,n))
- #define sqrmod(a,n) (((a) < HALF_WORD) ? ((a)*(a))%(n):_mulmod(a,a,n))
+ #define mulmod(a,b,n) _mulmod(a,b,n)
+ #define sqrmod(a,n) _mulmod(a,a,n)
#endif
diff --git a/primality.c b/primality.c
index 75dca8b..29804df 100644
--- a/primality.c
+++ b/primality.c
@@ -664,7 +664,7 @@ int _XS_is_almost_extra_strong_lucas_pseudoprime(UV n, UV increment)
int _XS_is_frobenius_underwood_pseudoprime(UV n)
{
int bit;
- UV x, result, a, b, np1, len, t1, t2, na;
+ UV x, result, a, b, np1, len, t1;
IV t;
if (n < 7) return (n == 2 || n == 3 || n == 5);
@@ -694,38 +694,26 @@ int _XS_is_frobenius_underwood_pseudoprime(UV n)
if (x == 0) {
result = mont5;
for (bit = len-2; bit >= 0; bit--) {
- t2 = addmod(b, b, n);
- na = mont_prod64(a, t2, n, npi);
- t1 = addmod(b, a, n);
- t2 = submod(b, a, n);
- b = mont_prod64(t1, t2, n, npi);
- a = na;
+ t1 = addmod(b, b, n);
+ b = mont_prod64(submod(b, a, n), addmod(b, a, n), n, npi);
+ a = mont_prod64(a, t1, n, npi);
if ( (np1 >> bit) & UVCONST(1) ) {
- t1 = addmod(a, a, n);
- na = addmod(t1, b, n);
- t1 = addmod(b, b, n);
- b = submod(t1, a, n);
- a = na;
+ t1 = b;
+ b = submod( addmod(b, b, n), a, n);
+ a = addmod( addmod(a, a, n), t1, n);
}
}
} else {
UV multiplier = addmod(x, mont2, n);
result = addmod( addmod(x, x, n), mont5, n);
for (bit = len-2; bit >= 0; bit--) {
- t1 = mont_prod64(a, x, n, npi);
- t2 = addmod(b, b, n);
- t1 = addmod(t1, t2, n);
- na = mont_prod64(a, t1, n, npi);
- t1 = addmod(b, a, n);
- t2 = submod(b, a, n);
- b = mont_prod64(t1, t2, n, npi);
- a = na;
+ t1 = addmod( mont_prod64(a, x, n, npi), addmod(b, b, n), n);
+ b = mont_prod64(submod(b, a, n), addmod(b, a, n), n, npi);
+ a = mont_prod64(a, t1, n, npi);
if ( (np1 >> bit) & UVCONST(1) ) {
- t1 = mont_prod64(a, multiplier, n, npi);
- na = addmod(t1, b, n);
- t1 = addmod(b, b, n);
- b = submod(t1, a, n);
- a = na;
+ t1 = b;
+ b = submod( addmod(b, b, n), a, n);
+ a = addmod( mont_prod64(a, multiplier, n, npi), t1, n);
}
}
}
@@ -738,38 +726,26 @@ int _XS_is_frobenius_underwood_pseudoprime(UV n)
if (x == 0) {
result = 5;
for (bit = len-2; bit >= 0; bit--) {
- t2 = addmod(b, b, n);
- na = mulmod(a, t2, n);
- t1 = addmod(b, a, n);
- t2 = submod(b, a, n);
- b = mulmod(t1, t2, n);
- a = na;
+ t1 = addmod(b, b, n);
+ b = mulmod( submod(b, a, n), addmod(b, a, n), n);
+ a = mulmod(a, t1, n);
if ( (np1 >> bit) & UVCONST(1) ) {
- t1 = addmod(a, a, n);
- na = addmod(t1, b, n);
- t1 = addmod(b, b, n);
- b = submod(t1, a, n);
- a = na;
+ t1 = b;
+ b = submod( addmod(b, b, n), a, n);
+ a = addmod( addmod(a, a, n), t1, n);
}
}
} else {
UV multiplier = addmod(x, 2, n);
result = addmod( addmod(x, x, n), 5, n);
for (bit = len-2; bit >= 0; bit--) {
- t1 = mulmod(a, x, n);
- t2 = addmod(b, b, n);
- t1 = addmod(t1, t2, n);
- na = mulmod(a, t1, n);
- t1 = addmod(b, a, n);
- t2 = submod(b, a, n);
- b = mulmod(t1, t2, n);
- a = na;
+ t1 = addmod( mulmod(a, x, n), addmod(b, b, n), n);
+ b = mulmod(submod(b, a, n), addmod(b, a, n), n);
+ a = mulmod(a, t1, n);
if ( (np1 >> bit) & UVCONST(1) ) {
- t1 = mulmod(a, multiplier, n);
- na = addmod(t1, b, n);
- t1 = addmod(b, b, n);
- b = submod(t1, a, n);
- a = na;
+ t1 = b;
+ b = submod( addmod(b, b, n), a, n);
+ a = addmod( mulmod(a, multiplier, n), t1, n);
}
}
}
@@ -798,7 +774,7 @@ static const UV mr_bases_large_3[3] = { UVCONST( 4230279247111683200 ),
static const UV mr_bases_large_7[7] = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 };
#endif
-int _XS_is_prob_prime(UV n)
+int is_prob_prime(UV n)
{
int ret;
diff --git a/primality.h b/primality.h
index ec7e0f9..1ac9744 100644
--- a/primality.h
+++ b/primality.h
@@ -17,6 +17,6 @@ extern int _XS_is_frobenius_underwood_pseudoprime(UV n);
extern int _XS_is_almost_extra_strong_lucas_pseudoprime(UV n, UV increment);
extern int _XS_BPSW(UV const n);
-extern int _XS_is_prob_prime(UV n);
+extern int is_prob_prime(UV n);
#endif
diff --git a/sieve.c b/sieve.c
index f359747..c28fc58 100644
--- a/sieve.c
+++ b/sieve.c
@@ -210,7 +210,7 @@ static void memtile(unsigned char* src, UV from, UV to) {
static UV sieve_prefill(unsigned char* mem, UV startd, UV endd)
{
- UV next_prime = 17;
+ UV vnext_prime = 17;
UV nbytes = endd - startd + 1;
MPUassert( (mem != 0) && (endd >= startd), "sieve_prefill bad arguments");
@@ -228,7 +228,7 @@ static UV sieve_prefill(unsigned char* mem, UV startd, UV endd)
if (startd == 0) mem[0] = 0x01; /* Correct first byte */
}
/* Leaving option open to tile 17 out and sieve, then return 19 */
- return next_prime;
+ return vnext_prime;
}
/* Wheel 30 sieve. Ideas from Terje Mathisen and Quesada / Van Pelt. */
diff --git a/util.c b/util.c
index ec16748..8480572 100644
--- a/util.c
+++ b/util.c
@@ -151,29 +151,6 @@ static UV count_zero_bits(const unsigned char* m, UV nbytes)
}
-/* Does trial division or prob tests, assuming x not divisible by 2, 3, or 5 */
-static int _is_prime7(UV n)
-{
- UV limit, i;
-
- if (n > MPU_PROB_PRIME_BEST)
- return _XS_is_prob_prime(n); /* We know this works for all 64-bit n */
-
- limit = isqrt(n);
- i = 7;
- while (1) { /* trial division, skipping multiples of 2/3/5 */
- if (i > limit) break; if ((n % i) == 0) return 0; i += 4;
- if (i > limit) break; if ((n % i) == 0) return 0; i += 2;
- if (i > limit) break; if ((n % i) == 0) return 0; i += 4;
- if (i > limit) break; if ((n % i) == 0) return 0; i += 2;
- if (i > limit) break; if ((n % i) == 0) return 0; i += 4;
- if (i > limit) break; if ((n % i) == 0) return 0; i += 6;
- if (i > limit) break; if ((n % i) == 0) return 0; i += 2;
- if (i > limit) break; if ((n % i) == 0) return 0; i += 6;
- }
- return 2;
-}
-
/* We'll use this little static sieve to quickly answer small values of
* is_prime, next_prime, prev_prime, prime_count
@@ -237,11 +214,11 @@ int _XS_is_prime(UV n)
if (isprime >= 0)
return isprime;
}
- return _XS_is_prob_prime(n);
+ return is_prob_prime(n);
}
-UV _XS_next_prime(UV n)
+UV next_prime(UV n)
{
UV d, m, sieve_size, next;
const unsigned char* sieve;
@@ -263,7 +240,7 @@ UV _XS_next_prime(UV n)
/* Move forward one, knowing we may not be on the wheel */
if (m == 29) { d++; m = 1; } else { m = nextwheel30[m]; }
n = d*30+m;
- while (!_is_prime7(n)) {
+ while (!is_prob_prime(n)) {
/* Move forward one, knowing we are on the wheel */
n += wheeladvance30[m];
m = nextwheel30[m];
@@ -272,7 +249,7 @@ UV _XS_next_prime(UV n)
}
-UV _XS_prev_prime(UV n)
+UV prev_prime(UV n)
{
const unsigned char* sieve;
UV d, m, prev;
@@ -293,7 +270,7 @@ UV _XS_prev_prime(UV n)
m = prevwheel30[m];
if (m==29) d--;
n = d*30+m;
- } while (!_is_prime7(n));
+ } while (!is_prob_prime(n));
return n;
}
@@ -704,7 +681,7 @@ UV _XS_nth_prime(UV n)
if (count >= n) { /* Too far. Walk backwards */
if (_XS_is_prime(lower_limit)) count--;
for (p = 0; p <= (count-n); p++)
- lower_limit = _XS_prev_prime(lower_limit);
+ lower_limit = prev_prime(lower_limit);
return lower_limit;
}
count -= 3;
@@ -753,12 +730,10 @@ signed char* _moebius_range(UV lo, UV hi)
/* Kuznetsov indicates that the Deléglise & Rivat (1996) method can be
* modified to work on logs, which allows us to operate with no
* intermediate memory at all. Same time as the D&R method, less memory. */
- unsigned char* A;
unsigned char logp;
UV nextlog;
Newz(0, mu, hi-lo+1, signed char);
- A = (unsigned char*) mu;
if (sqrtn*sqrtn != hi) sqrtn++; /* ceil sqrtn */
logp = 1; nextlog = 3; /* 2+1 */
@@ -769,19 +744,20 @@ signed char* _moebius_range(UV lo, UV hi)
nextlog = ((nextlog-1)*4)+1;
}
for (i = PGTLO(p, lo); i <= hi; i += p)
- A[i-lo] += logp;
+ mu[i-lo] += logp;
for (i = PGTLO(p2, lo); i <= hi; i += p2)
- A[i-lo] |= 0x80;
+ mu[i-lo] |= 0x80;
} END_DO_FOR_EACH_PRIME
logp = log2floor(lo);
nextlog = 2UL << logp;
for (i = lo; i <= hi; i++) {
- unsigned char a = A[i-lo];
+ unsigned char a = mu[i-lo];
if (i >= nextlog) { logp++; nextlog *= 2; } /* logp is log(p)/log(2) */
- if (a & 0x80) { mu[i-lo] = 0; }
- else if (a >= logp) { mu[i-lo] = 1 - 2*(a&1); }
- else { mu[i-lo] = -1 + 2*(a&1); }
+ if (a & 0x80) { a = 0; }
+ else if (a >= logp) { a = 1 - 2*(a&1); }
+ else { a = -1 + 2*(a&1); }
+ mu[i-lo] = a;
}
if (lo == 0) mu[0] = 0;
@@ -1006,15 +982,13 @@ UV znorder(UV a, UV n) {
phi = carmichael_lambda(n);
nfactors = factor_exp(phi, fac, exp);
for (i = 0; i < nfactors; i++) {
- UV b, ek;
- UV pi = fac[i], ei = exp[i];
+ UV b, ek, pi = fac[i], ei = exp[i];
UV phidiv = phi / pi;
- for (j = 1; j < ei; j++) phidiv /= pi;
+ for (j = 1; j < ei; j++)
+ phidiv /= pi;
b = powmod(a, phidiv, n);
- ek = 0;
- while (b != 1) {
+ for (ek = 0; b != 1; b = powmod(b, pi, n)) {
if (ek++ >= ei) return 0;
- b = powmod(b, pi, n);
k *= pi;
}
}
@@ -1030,7 +1004,7 @@ UV znprimroot(UV n) {
if (n % 4 == 0) return 0;
phi = totient(n);
/* Check if a primitive root exists. */
- if (!_XS_is_prob_prime(n) && phi != carmichael_lambda(n)) return 0;
+ if (!is_prob_prime(n) && phi != carmichael_lambda(n)) return 0;
nfactors = factor_exp(phi, fac, exp);
for (i = 0; i < nfactors; i++)
exp[i] = phi / fac[i]; /* exp[i] = phi(n) / i-th-factor-of-phi(n) */
diff --git a/util.h b/util.h
index ed3c8f2..e92e32c 100644
--- a/util.h
+++ b/util.h
@@ -9,8 +9,8 @@ extern int _XS_get_callgmp(void);
extern void _XS_set_callgmp(int v);
extern int _XS_is_prime(UV x);
-extern UV _XS_next_prime(UV x);
-extern UV _XS_prev_prime(UV x);
+extern UV next_prime(UV x);
+extern UV prev_prime(UV x);
extern UV _XS_prime_count(UV low, UV high);
extern UV _XS_nth_prime(UV x);
@@ -41,14 +41,6 @@ extern UV znprimroot(UV n);
extern UV znorder(UV a, UV n);
extern UV znlog(UV a, UV g, UV p);
-/* Above this value, is_prime will do deterministic Miller-Rabin */
-/* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */
-#if (BITS_PER_WORD == 64) || HAVE_STD_U64
- #define MPU_PROB_PRIME_BEST UVCONST(100000)
-#else
- #define MPU_PROB_PRIME_BEST UVCONST(100000000)
-#endif
-
#if defined(FUNC_isqrt)
static UV isqrt(UV n) {
UV root;
--
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