[libmath-prime-util-perl] 15/20: Move totient range to util.c, speed it up a little bit

[Partha P. Mukherjee]

This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.23
in repository libmath-prime-util-perl.

commit 4af68950730378cac15f16f221289ec81146b70e
Author: Dana Jacobsen <dana at acm.org>
Date:   Mon Mar 4 05:25:35 2013 -0800

    Move totient range to util.c, speed it up a little bit
---
 XS.xs  | 39 ++++++++++++---------------------------
 util.c | 32 +++++++++++++++++++++++++++++++-
 util.h |  1 +
 3 files changed, 44 insertions(+), 28 deletions(-)

diff --git a/XS.xs b/XS.xs
index bd60751..af95cde 100644
--- a/XS.xs
+++ b/XS.xs
@@ -359,46 +359,31 @@ _XS_RiemannR(double x)
 
 void
 _XS_totient(IN UV lo, IN UV hi = 0)
+  PREINIT:
+    UV i;
   PPCODE:
-    if (hi != lo && hi != 0) {
-      /* Totients in a range, returns array */
+    if (hi != lo && hi != 0) { /* Totients in a range, returns array */
       UV* totients;
-      UV i, p;
-
       if (hi < lo) XSRETURN_EMPTY;
       if (lo < 2) {
         if (lo <= 0           ) XPUSHs(sv_2mortal(newSVuv(0)));
         if (lo <= 1 && hi >= 1) XPUSHs(sv_2mortal(newSVuv(1)));
         lo = 2;
       }
-      New(0, totients, hi-lo+1, UV);
-      if (totients == 0)
-        croak("Could not get memory for %"UVuf" totients\n", hi);
-      for (i = lo; i <= hi; i++)
-        totients[i-lo] = i;
-      prime_precalc( hi/2 );
-      for (p = 2; p <= hi/2; p = _XS_next_prime(p)) {
-        i = 2*p;
-        if (i < lo)  i = p*(lo/p) + ( (lo%p) ? p : 0 );
-        for ( ; i <= hi; i += p)
-          totients[i-lo] -= totients[i-lo]/p;
+      if (hi >= lo) {
+        totients = _totient_range(lo, hi);
+        /* Extend the stack to handle how many items we'll return */
+        EXTEND(SP, hi-lo+1);
+        for (i = lo; i <= hi; i++)
+          PUSHs(sv_2mortal(newSVuv(totients[i-lo])));
+        Safefree(totients);
       }
-      /* Extend the stack to handle how many items we'll return */
-      EXTEND(SP, hi-lo+1);
-      for (i = lo; i <= hi; i++) {
-        UV t = totients[i-lo];
-        if (t == i)
-          t = i-1;
-        PUSHs(sv_2mortal(newSVuv(t)));
-      }
-      Safefree(totients);
-
     } else {
       UV facs[MPU_MAX_FACTORS+1];  /* maximum number of factors is log2n */
-      UV i, nfacs, totient, lastf;
+      UV nfacs, totient, lastf;
       UV n = lo;
       if (n <= 1) XSRETURN_UV(n);
-      nfacs = trial_factor(n, facs, 0);
+      nfacs = factor(n, facs);
       totient = 1;
       lastf = 0;
       for (i = 0; i < nfacs; i++) {
diff --git a/util.c b/util.c
index 66525c1..1da3426 100644
--- a/util.c
+++ b/util.c
@@ -651,7 +651,7 @@ UV _XS_nth_prime(UV n)
 
 /* Return an IV array with lo-hi+1 elements.  mu[k-lo] = µ(k) for k = lo .. hi.
  * It is the callers responsibility to call Safefree on the result. */
-#define PGTLO(p,lo)  (p >= lo) ? p : (p*(lo/p) + ((lo%p)?p:0))
+#define PGTLO(p,lo)  ((p) >= lo) ? (p) : ((p)*(lo/(p)) + ((lo%(p))?(p):0))
 char* _moebius_range(UV lo, UV hi)
 {
   char* mu;
@@ -746,6 +746,36 @@ char* _moebius_range(UV lo, UV hi)
   return mu;
 }
 
+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);
+  if (totients == 0)
+    croak("Could not get memory for %"UVuf" totients\n", hi);
+  for (i = lo; i <= hi; i++)
+    totients[i-lo] = i;
+  for (i = PGTLO(2*2, lo); i <= hi; i += 2) totients[i-lo] -= totients[i-lo]/2;
+  for (i = PGTLO(2*3, lo); i <= hi; i += 3) totients[i-lo] -= totients[i-lo]/3;
+  for (i = PGTLO(2*5, lo); i <= hi; i += 5) totients[i-lo] -= totients[i-lo]/5;
+  sievehi = hi/2;
+  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 = PGTLO(2*p, lo); i <= hi; i += p)
+        totients[i-lo] -= totients[i-lo]/p;
+    } END_DO_FOR_EACH_SIEVE_PRIME
+    release_prime_cache(sieve);
+  }
+  for (i = lo; i <= hi; i++)
+    if (totients[i-lo] == i)
+      totients[i-lo] = i-1;
+  return totients;
+}
+
 IV _XS_mertens(UV n) {
 #if 0
   /* Benito and Varona 2008, theorem 3.  Segment. */
diff --git a/util.h b/util.h
index a7b56e4..a2d7f7c 100644
--- a/util.h
+++ b/util.h
@@ -15,6 +15,7 @@ extern UV  _XS_prime_count(UV low, UV high);
 extern UV  _XS_nth_prime(UV x);
 
 extern char*  _moebius_range(UV low, UV high);
+extern UV*    _totient_range(UV low, UV high);
 extern IV     _XS_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