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 517233e310bf4209ce5b25176941c0cd19be715f
Author: Dana Jacobsen <d...@acm.org>
Date:   Tue Dec 31 01:31:24 2013 -0800

    More XS changes, plus small optimizations from bulk88
---
 XS.xs                  | 117 ++++++++++++++++++++-----------------------------
 aks.c                  |   6 ++-
 lib/Math/Prime/Util.pm |   8 ++--
 util.c                 |  35 +++++++--------
 util.h                 |   2 +-
 5 files changed, 73 insertions(+), 95 deletions(-)

diff --git a/XS.xs b/XS.xs
index 851142f..2d8e507 100644
--- a/XS.xs
+++ b/XS.xs
@@ -129,9 +129,11 @@ static int _validate_int(pTHX_ SV* n, int negok)
 static int _vcallsubn(pTHX_ I32 flags, const char* name, int nargs)
 {
     dSP;
+    char fullname[80] = "Math::Prime::Util::";
+    strncat(fullname, name, 60);
     PUSHMARK(SP-nargs);
     PUTBACK;
-    return call_pv(name, flags);
+    return call_pv(fullname, flags);
 }
 #define _vcallsub(func) (void)_vcallsubn(aTHX_ G_SCALAR, func, 1)
 
@@ -163,7 +165,7 @@ prime_memfree()
       case 2:  ret = _XS_get_verbose(); break;
       case 3:  ret = _XS_get_callgmp(); break;
       case 4:  ret = get_prime_cache(0,0); break;
-      case 5:  ret = BITS_PER_WORD;     break;
+      default: ret = BITS_PER_WORD;     break;
     }
     if (ix > 1) XSRETURN_UV(ret);
 
@@ -223,7 +225,6 @@ _XS_nth_prime(IN UV n)
 UV
 _XS_legendre_phi(IN UV x, IN UV a)
 
-
 SV*
 sieve_primes(IN UV low, IN UV high)
   ALIAS:
@@ -356,9 +357,6 @@ trial_factor(IN UV n, ...)
     }
 
 int
-_XS_is_pseudoprime(IN UV n, IN UV a)
-
-int
 _XS_miller_rabin(IN UV n, ...)
   PREINIT:
     UV bases[64];
@@ -396,15 +394,17 @@ _XS_lucas_sequence(IN UV n, IN IV P, IN IV Q, IN UV k)
     XPUSHs(sv_2mortal(newSVuv( Qk )));
 
 int
-_XS_is_prime(IN UV n)
+_XS_is_prime(IN UV n, ...)
   ALIAS:
     _XS_is_prob_prime = 1
     _XS_is_lucas_pseudoprime = 2
     _XS_is_strong_lucas_pseudoprime = 3
     _XS_is_extra_strong_lucas_pseudoprime = 4
     _XS_is_frobenius_underwood_pseudoprime = 5
-    _XS_is_aks_prime = 6
-    _XS_BPSW = 7
+    _XS_is_almost_extra_strong_lucas_pseudoprime = 6
+    _XS_is_pseudoprime = 7
+    _XS_is_aks_prime = 8
+    _XS_BPSW = 9
   PREINIT:
     int ret;
   CODE:
@@ -415,19 +415,16 @@ _XS_is_prime(IN UV n)
       case 3: ret = _XS_is_lucas_pseudoprime(n, 1); break;
       case 4: ret = _XS_is_lucas_pseudoprime(n, 2); break;
       case 5: ret = _XS_is_frobenius_underwood_pseudoprime(n); break;
-      case 6: ret = _XS_is_aks_prime(n); break;
+      case 6: ret = _XS_is_almost_extra_strong_lucas_pseudoprime
+                    ( n, (items == 1) ? 1 : my_svuv(ST(1)) );  break;
+      case 7: ret = _XS_is_pseudoprime(n, my_svuv(ST(1)));     break;
+      case 8: ret = _XS_is_aks_prime(n); break;
       default:ret = _XS_BPSW(n); break;
     }
     RETVAL = ret;
   OUTPUT:
     RETVAL
 
-int
-_XS_is_almost_extra_strong_lucas_pseudoprime(IN UV n, IN UV increment = 1)
-  CODE:
-    RETVAL = _XS_is_almost_extra_strong_lucas_pseudoprime(n, increment);
-  OUTPUT:
-    RETVAL
 
 void
 is_prime(IN SV* svn)
@@ -445,11 +442,11 @@ is_prime(IN SV* svn)
     } else {
       const char* sub = 0;
       if (_XS_get_callgmp())
-        sub = (ix == 0) ? "Math::Prime::Util::GMP::is_prime"
-                        : "Math::Prime::Util::GMP::is_prob_prime";
+        sub = (ix == 0) ? "GMP::is_prime"
+                        : "GMP::is_prob_prime";
       else
-        sub = (ix == 0) ? "Math::Prime::Util::_generic_is_prime"
-                        : "Math::Prime::Util::_generic_is_prob_prime";
+        sub = (ix == 0) ? "_generic_is_prime"
+                        : "_generic_is_prob_prime";
       _vcallsub(sub);
       return; /* skip implicit PUTBACK */
     }
@@ -468,8 +465,8 @@ next_prime(IN SV* svn)
           XSRETURN_UV(_XS_next_prime(n));
       }
     }
-    _vcallsub((ix == 0) ?  "Math::Prime::Util::_generic_next_prime" :
-                           "Math::Prime::Util::_generic_prev_prime" );
+    _vcallsub((ix == 0) ?  "_generic_next_prime" :
+                           "_generic_prev_prime" );
     return; /* skip implicit PUTBACK */
 
 void
@@ -489,7 +486,7 @@ factor(IN SV* svn)
         }
       }
     } else {
-      _vcallsubn(aTHX_ GIMME_V, "Math::Prime::Util::_generic_factor", 1);
+      _vcallsubn(aTHX_ GIMME_V, "_generic_factor", 1);
       return; /* skip implicit PUTBACK */
     }
 
@@ -505,7 +502,7 @@ divisor_sum(IN SV* svn, ...)
       UV sigma = divisor_sum(n, k);
       if (sigma != 0)  XSRETURN_UV(sigma);   /* sigma 0 means overflow */
     }
-    _vcallsubn(aTHX_ G_SCALAR, 
"Math::Prime::Util::_generic_divisor_sum",items);
+    _vcallsubn(aTHX_ G_SCALAR, "_generic_divisor_sum",items);
     return; /* skip implicit PUTBACK */
 
 
@@ -521,7 +518,7 @@ znorder(IN SV* sva, IN SV* svn)
       if (order == 0) XSRETURN_UNDEF;
       XSRETURN_UV(order);
     }
-    _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_znorder", 2);
+    _vcallsubn(aTHX_ G_SCALAR, "_generic_znorder", 2);
     return; /* skip implicit PUTBACK */
 
 void
@@ -536,7 +533,7 @@ znprimroot(IN SV* svn)
       if (r == 0 && n != 1) XSRETURN_EMPTY;
       XSRETURN_UV(r);
     }
-    _vcallsub("Math::Prime::Util::_generic_znprimroot");
+    _vcallsub("_generic_znprimroot");
     return; /* skip implicit PUTBACK */
 
 void
@@ -554,7 +551,7 @@ kronecker(IN SV* sva, IN SV* svb)
       IV b = my_sviv(svb);
       XSRETURN_IV( kronecker_ss(a, b) );
     }
-    _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_kronecker", 2);
+    _vcallsubn(aTHX_ G_SCALAR, "_generic_kronecker", 2);
     return; /* skip implicit PUTBACK */
 
 double
@@ -591,7 +588,7 @@ euler_phi(IN SV* svlo, ...)
         UV lo = my_svuv(svlo);
         XSRETURN_UV(totient(lo));
       } else {
-        _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_euler_phi", 1);
+        _vcallsubn(aTHX_ G_SCALAR, "_generic_euler_phi", 1);
         return; /* skip implicit PUTBACK */
       }
     } else if (items == 2) {
@@ -617,7 +614,7 @@ euler_phi(IN SV* svlo, ...)
           Safefree(totients);
         }
       } else {
-        _vcallsubn(aTHX_ 
G_ARRAY,"Math::Prime::Util::_generic_euler_phi",items);
+        _vcallsubn(aTHX_ G_ARRAY,"_generic_euler_phi",items);
         return; /* skip implicit PUTBACK */
       }
     } else {
@@ -633,7 +630,7 @@ moebius(IN SV* svlo, ...)
         UV n = my_svuv(svlo);
         XSRETURN_IV(moebius(n));
       } else {
-        _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_moebius",1);
+        _vcallsubn(aTHX_ G_SCALAR, "_generic_moebius",1);
         return; /* skip implicit PUTBACK */
       }
     } else if (items == 2) {
@@ -655,34 +652,32 @@ moebius(IN SV* svlo, ...)
           Safefree(mu);
         }
       } else {
-        _vcallsubn(aTHX_ G_ARRAY, "Math::Prime::Util::_generic_moebius",items);
+        _vcallsubn(aTHX_ G_ARRAY, "_generic_moebius",items);
         return; /* skip implicit PUTBACK */
       }
     } else {
       croak("Usage: moebius(n) or moebius(1o,hi)");
     }
 
-IV
-_XS_mertens(IN UV n)
-
-
 void
 carmichael_lambda(IN SV* svn)
   ALIAS:
-    exp_mangoldt = 1
-  PREINIT:
-    int status;
+    mertens = 1
+    exp_mangoldt = 2
   PPCODE:
-    status = _validate_int(aTHX_ svn, (ix == 0) ? 0 : 1);
-    if (status == -1) {
-      XSRETURN_UV(1);
-    } else if (status == 1) {
-      UV n = my_svuv(svn);
-      if (ix == 0) XSRETURN_UV(carmichael_lambda(n));
-      else         XSRETURN_UV(exp_mangoldt(n));
+    int status = _validate_int(aTHX_ svn, (ix > 1) ? 1 : 0);
+    switch (ix) {
+      case 0: if (status == 1) XSRETURN_UV(carmichael_lambda(my_svuv(svn)));
+              _vcallsub("_generic_carmichael_lambda");
+              break;
+      case 1: if (status == 1) XSRETURN_IV(mertens(my_svuv(svn)));
+              _vcallsub("_generic_mertens");
+              break;
+      default:if (status ==-1) XSRETURN_UV(1);
+              if (status == 1) XSRETURN_UV(exp_mangoldt(my_svuv(svn)));
+              _vcallsub("_generic_exp_mangoldt");
+              break;
     }
-    _vcallsub( (ix == 0) ? "Math::Prime::Util::_generic_carmichael_lambda"
-                         : "Math::Prime::Util::_generic_exp_mangoldt" );
     return; /* skip implicit PUTBACK */
 
 int
@@ -717,11 +712,7 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
   PPCODE:
   {
 #if !USE_MULTICALL
-    PUSHMARK(SP);
-    XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
-    PUTBACK;
-    (void) call_pv("Math::Prime::Util::_generic_forprimes", G_VOID|G_DISCARD);
-    SPAGAIN;
+    _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_forprimes",items);
 #else
     UV beg, end;
     GV *gv;
@@ -733,11 +724,7 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
     UV seg_base, seg_low, seg_high;
 
     if (!_validate_int(aTHX_ svbeg, 0) || (items >= 3 && !_validate_int(aTHX_ 
svend,0))) {
-      PUSHMARK(SP);
-      XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
-      PUTBACK;
-      (void) call_pv("Math::Prime::Util::_generic_forprimes", 
G_VOID|G_DISCARD);
-      SPAGAIN;
+      _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_forprimes",items);
       return;
     }
     if (items < 3) {
@@ -828,11 +815,7 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
     if (items <= 1) return;
 
     if (!_validate_int(aTHX_ svbeg, 0) || (items >= 3 && !_validate_int(aTHX_ 
svend,0))) {
-      PUSHMARK(SP);
-      XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
-      PUTBACK;
-      (void) call_pv("Math::Prime::Util::_generic_forcomposites", 
G_VOID|G_DISCARD);
-      SPAGAIN;
+      _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_forcomposites",items);
       return;
     }
 
@@ -868,9 +851,9 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
       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 ) {
-          cbeg = (prevprime+1 < beg)   ? beg : prevprime+1;
+          cbeg = prevprime+1;  if (cbeg < beg) cbeg = beg;
           prevprime = seg_base + p;
-          cend = (prevprime-1 > end)   ? end : prevprime-1;
+          cend = prevprime-1;  if (cend > end) cbeg = end;
           for (c = cbeg; c <= cend; c++) {
             sv_setuv(svarg, c);  MULTICALL;
           }
@@ -913,11 +896,7 @@ fordivisors (SV* block, IN SV* svn)
     if (items <= 1) return;
 
     if (!_validate_int(aTHX_ svn, 0)) {
-      PUSHMARK(SP);
-      XPUSHs(block); XPUSHs(svn);
-      PUTBACK;
-      (void) call_pv("Math::Prime::Util::_generic_fordivisors", 
G_VOID|G_DISCARD);
-      SPAGAIN;
+      _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_fordivisors",2);
       return;
     }
 
diff --git a/aks.c b/aks.c
index 6f4f9d8..733cb11 100644
--- a/aks.c
+++ b/aks.c
@@ -43,7 +43,6 @@ static int is_perfect_power(UV n) {
   UV b, last;
   if ((n <= 3) || (n == UV_MAX)) return 0;
   if ((n & (n-1)) == 0)  return 1;          /* powers of 2    */
-  last = log2floor(n-1) + 1;
 #if (BITS_PER_WORD == 32) || (DBL_DIG > 19)
   if (1) {
 #elif DBL_DIG == 10
@@ -54,13 +53,16 @@ static int is_perfect_power(UV n) {
   if ( n < (UV) pow(10, DBL_DIG) ) {
 #endif
     /* Simple floating point method.  Fast, but need enough mantissa. */
-    b = isqrt(n); if (b*b == n)  return 1; /* perfect square */
+    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)) {
       UV root = (UV) (pow(n, 1.0 / (double)b) + 0.5);
       if ( ((UV)(pow(root, b)+0.5)) == n)  return 1;
     }
   } else {
     /* Dietzfelbinger, algorithm 2.3.5 (without optimized exponential) */
+    last = log2floor(n-1) + 1;
     for (b = 2; b <= last; b++) {
       UV a = 1;
       UV c = n;
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 28f5a19..44c5b4d 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -66,8 +66,6 @@ sub _import_nobigint {
   undef *nth_prime;       *nth_prime         = \&_XS_nth_prime;
   undef *is_pseudoprime;  *is_pseudoprime    = \&_XS_is_pseudoprime;
   undef *is_strong_pseudoprime;  *is_strong_pseudoprime = \&_XS_miller_rabin;
-  undef *moebius;         *moebius           = \&_XS_moebius;
-  undef *mertens;         *mertens           = \&_XS_mertens;
   undef *chebyshev_theta; *chebyshev_theta   = \&_XS_chebyshev_theta;
   undef *chebyshev_psi;   *chebyshev_psi     = \&_XS_chebyshev_psi;
   # These should be fast anyway, but this skips validation.
@@ -105,6 +103,7 @@ BEGIN {
     *exp_mangoldt  = \&Math::Prime::Util::_generic_exp_mangoldt;
     *euler_phi     = \&Math::Prime::Util::_generic_euler_phi;
     *moebius       = \&Math::Prime::Util::_generic_moebius;
+    *mertens       = \&Math::Prime::Util::_generic_mertens;
     *carmichael_lambda = \&Math::Prime::Util::_generic_carmichael_lambda;
     *kronecker     = \&Math::Prime::Util::_generic_kronecker;
     *divisor_sum   = \&Math::Prime::Util::_generic_divisor_sum;
@@ -1275,10 +1274,9 @@ sub _generic_moebius {
 }
 
 # A002321 Mertens' function.  mertens(n) = sum(moebius(1,n))
-sub mertens {
+sub _generic_mertens {
   my($n) = @_;
   _validate_num($n) || _validate_positive_integer($n);
-  return _XS_mertens($n) if $n <= $_XS_MAXVAL;
   # This is the most basic Del├ęglise and Rivat algorithm.  u = n^1/2
   # and no segmenting is done.  Their algorithm uses u = n^1/3, breaks
   # the summation into two parts, and calculates those in segments.  Their
@@ -4343,7 +4341,7 @@ Project Euler, problem 21 (Amicable numbers):
 
 Project Euler, problem 41 (Pandigital prime), brute force command line:
 
-  perl -MMath::Prime::Util=:all -E 'my @p = grep { 
/1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} @{primes(1000000,9999999)}; say $p[-1];'
+  perl -MMath::Prime::Util=primes -MList::Util=first -E 'say first { 
/1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} reverse @{primes(1000000,9999999)};'
 
 Project Euler, problem 47 (Distinct primes factors):
 
diff --git a/util.c b/util.c
index 4f14cae..78d1714 100644
--- a/util.c
+++ b/util.c
@@ -197,17 +197,15 @@ static const unsigned char prime_sieve30[] =
 /* Return of 2 if n is prime, 0 if not.  Do it fast. */
 int _XS_is_prime(UV n)
 {
+  if (n <= 10)
+    return (n == 2 || n == 3 || n == 5 || n == 7) ? 2 : 0;
+
   if (n < UVCONST(2000000000)) {
-    UV d, m;
-    unsigned char mtab;
+    UV d = n/30;
+    UV m = n - d*30;
+    unsigned char mtab = masktab30[m];  /* Bitmask in mod30 wheel */
     const unsigned char* sieve;
-    int isprime = -1;
-
-    if (n <= 10)  return (n == 2 || n == 3 || n == 5 || n == 7) ? 2 : 0;
-
-    d = n/30;
-    m = n - d*30;
-    mtab = masktab30[m];  /* Bitmask in mod30 wheel */
+    int isprime;
 
     /* Return 0 if a multiple of 2, 3, or 5 */
     if (mtab == 0)
@@ -220,6 +218,7 @@ int _XS_is_prime(UV n)
     if (!(n%7) || !(n%11) || !(n%13)) return 0;
 
     /* Check primary cache */
+    isprime = -1;
     if (n <= get_prime_cache(0, &sieve))
       isprime = 2*((sieve[d] & mtab) == 0);
     release_prime_cache(sieve);
@@ -635,12 +634,12 @@ static const unsigned short primes_small[] =
 /* The nth prime will be less or equal to this number */
 static UV _XS_nth_prime_upper(UV n)
 {
-  double fn = (double) n;
-  double flogn, flog2n, upper;
+  double fn, flogn, flog2n, upper;
 
   if (n < NPRIMES_SMALL)
     return primes_small[n];
 
+  fn     = (double) n;
   flogn  = log(n);
   flog2n = log(flogn);    /* Note distinction between log_2(n) and log^2(n) */
 
@@ -846,7 +845,7 @@ UV* _totient_range(UV lo, UV hi) {
   return totients;
 }
 
-IV _XS_mertens(UV n) {
+IV mertens(UV n) {
   /* See Del├ęglise and Rivat (1996) for O(n^2/3 log(log(n))^1/3) algorithm.
    * This implementation uses their lemma 2.1 directly, so is ~ O(n).
    * In serial it is quite a bit faster than segmented summation of mu
@@ -948,11 +947,11 @@ int kronecker_ss(IV a, IV b) {
 }
 
 UV totient(UV n) {
-  UV i, facs[MPU_MAX_FACTORS+1];
-  UV nfacs = factor(n, facs);
-  UV totient = 1;
-  UV lastf = 0;
+  UV i, nfacs, totient, lastf, facs[MPU_MAX_FACTORS+1];
   if (n <= 1) return n;
+  nfacs = factor(n, facs);
+  totient = 1;
+  lastf = 0;
   for (i = 0; i < nfacs; i++) {
     UV f = facs[i];
     if (f == lastf) { totient *= f;               }
@@ -1235,7 +1234,7 @@ double _XS_LogarithmicIntegral(double x) {
   if (x == 0) return 0;
   if (x == 1) return -INFINITY;
   if (x == 2) return li2;
-  if (x <= 0) croak("Invalid input to LogarithmicIntegral:  x must be > 0");
+  if (x < 0) croak("Invalid input to LogarithmicIntegral:  x must be >= 0");
   return _XS_ExponentialIntegral(log(x));
 }
 
@@ -1246,7 +1245,7 @@ UV _XS_Inverse_Li(UV x) {
   UV lo = (UV) (n*logn);
   UV hi = (UV) (n*logn * 2 + 2);
 
-  if (x < 1)  return 0;
+  if (x == 0)  return 0;
   if (hi <= lo) hi = UV_MAX;
   while (lo < hi) {
     UV mid = lo + (hi-lo)/2;
diff --git a/util.h b/util.h
index 6cb0cbb..0c40e59 100644
--- a/util.h
+++ b/util.h
@@ -17,7 +17,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     _XS_mertens(UV n);
+extern IV     mertens(UV n);
 extern double _XS_chebyshev_theta(UV n);
 extern double _XS_chebyshev_psi(UV n);
 

-- 
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