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 <d...@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

_______________________________________________
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

Reply via email to