[libmath-prime-util-perl] 40/43: More optimizations for euler_phi(0, n)

[Partha P. Mukherjee]

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

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

commit 7f6602eed22e6366958ebb4dfd094422fc826cec
Author: Dana Jacobsen <dana at acm.org>
Date:   Mon Apr 21 11:58:02 2014 -0700

    More optimizations for euler_phi(0,n)
---
 Changes |  2 +-
 TODO    |  3 ++-
 util.c  | 39 +++++++++++++++++++++++----------------
 3 files changed, 26 insertions(+), 18 deletions(-)

diff --git a/Changes b/Changes
index 3b12038..09ee7da 100644
--- a/Changes
+++ b/Changes
@@ -1,6 +1,6 @@
 Revision history for Perl module Math::Prime::Util
 
-0.40  2014-03
+0.40  2014-04-21
 
     [ADDED]
 
diff --git a/TODO b/TODO
index ae9e8e2..f8132c7 100644
--- a/TODO
+++ b/TODO
@@ -70,4 +70,5 @@
 - znlog better implementation
 
 - We don't use legendre_phi for other functions any more, but it'd be nice
-  to speed it up using some ideas from the Ohana 2011 SAGE code.
+  to speed it up using some ideas from the Ohana 2011 SAGE branch.  For example
+  (10**13,10**5) takes 2.5x longer, albeit with 6x less memory.
diff --git a/util.c b/util.c
index f9e8879..446e2fa 100644
--- a/util.c
+++ b/util.c
@@ -1071,10 +1071,9 @@ UV* _totient_range(UV lo, UV hi) {
     return totients;
   }
 
-  /* If doing a full sieve, do it monolithic.  More aux memory, but faster. */
+  /* If doing a full sieve, do it monolithic.  Faster. */
   if (lo == 0) {
     UV* prime;
-    char* setb;
     double loghi = log(hi);
     UV max_index = (hi < 67)     ? 18
                  : (hi < 355991) ? 15+(hi/(loghi-1.09))
@@ -1082,24 +1081,32 @@ UV* _totient_range(UV lo, UV hi) {
     UV j, index, nprimes = 0;
 
     New(0, prime, max_index, UV);  /* could use prime_count_upper(hi) */
-    Newz(0, setb, (hi+1+7)/8, char);
-    for (i = 2; i <= hi; i++) {
-      if ( (setb[i/8] & (1 << i%8)) == 0 ) {
-        totients[i] = i-1;
-        prime[nprimes++] = i;
-      }
-      for (j=0, index=2*i; j < nprimes && index <= hi; index = i*prime[++j]) {
-        setb[index/8] |= 1 << (index%8);
-        if (i % prime[j] == 0) {
-          totients[index] = totients[i]*prime[j];
-          break;
-        } else {
-          totients[index] = totients[i]*(prime[j]-1);
+    memset(totients, 0, (hi-lo+1) * sizeof(UV));
+    for (i = 2; i <= hi/2; i++) {
+      index = 2*i;
+      if ( !(i&1) ) {
+        if (i == 2) { totients[2] = 1; prime[nprimes++] = 2; }
+        totients[index] = totients[i]*2;
+      } else {
+        if (totients[i] == 0) {
+          totients[i] = i-1;
+          prime[nprimes++] = i;
+        }
+        for (j=0; j < nprimes && index <= hi; index = i*prime[++j]) {
+          if (i % prime[j] == 0) {
+            totients[index] = totients[i]*prime[j];
+            break;
+          } else {
+            totients[index] = totients[i]*(prime[j]-1);
+          }
         }
       }
     }
-    Safefree(setb);
     Safefree(prime);
+    /* All totient values have been filled in except the primes.  Mark them. */
+    for (i = ((hi/2) + 1) | 1; i <= hi; i += 2)
+      if (totients[i] == 0)
+        totients[i] = i-1;
     totients[1] = 1;
     totients[0] = 0;
     return totients;

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git