[libmath-prime-util-perl] 38/43: Add slightly faster ranged totient for start = 0

[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 6e9b8a26a2afc03bfc880846288cd5543b5c62e2
Author: Dana Jacobsen <dana at acm.org>
Date:   Fri Apr 18 15:16:29 2014 -0700

    Add slightly faster ranged totient for start = 0
---
 Changes             |  2 ++
 util.c              | 37 ++++++++++++++++++++++++++++++++++++-
 xt/totient-range.pl |  5 +++--
 3 files changed, 41 insertions(+), 3 deletions(-)

diff --git a/Changes b/Changes
index 45c1912..96765af 100644
--- a/Changes
+++ b/Changes
@@ -26,6 +26,8 @@ Revision history for Perl module Math::Prime::Util
     - Fix legendre_phi when given tiny x and large a (odd test case).
       Some speedups for huge a, from R. Andrew Ohana.
 
+    - Ranged totient is slightly faster with start of 0.
+
 
 0.39  2014-03-01
 
diff --git a/util.c b/util.c
index f246c34..f9e8879 100644
--- a/util.c
+++ b/util.c
@@ -583,6 +583,7 @@ UV twin_prime_count_approx(UV n)
   /* Best would be another estimate for n < ~ 5000 */
   if (n < 2000) return twin_prime_count(3,n);
   {
+    /* Sebah and Gourdon 2002 */
     const long double two_C2 = 1.32032363169373914785562422L;
     const long double two_over_log_two = 2.8853900817779268147198494L;
     long double ln = (long double) n;
@@ -1064,12 +1065,46 @@ UV* _totient_range(UV lo, UV hi) {
   New(0, totients, hi-lo+1, UV);
 
   /* Do via factoring if very small or if we have a small range */
-  if (hi < 100 || hi/(hi-lo+1) > 1000) {
+  if (hi < 100 || (hi-lo) < 10 || hi/(hi-lo+1) > 1000) {
     for (i = lo; i <= hi; i++)
       totients[i-lo] = totient(i);
     return totients;
   }
 
+  /* If doing a full sieve, do it monolithic.  More aux memory, but 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))
+                 : (hi/loghi) * (1.0+1.0/loghi+2.51/(loghi*loghi));
+    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);
+        }
+      }
+    }
+    Safefree(setb);
+    Safefree(prime);
+    totients[1] = 1;
+    totients[0] = 0;
+    return totients;
+  }
+
   for (i = lo; i <= hi; i++) {
     UV v = i;
     if (i % 2 == 0)  v -= v/2;
diff --git a/xt/totient-range.pl b/xt/totient-range.pl
index 6371a6c..035973b 100755
--- a/xt/totient-range.pl
+++ b/xt/totient-range.pl
@@ -14,9 +14,10 @@ print "...";
 unshift @phi, 0;
 print "...done\n";
 
+print "Running non-stop random tests.  Break when desired.\n";
 while (1) {
-  my $beg = 1 + int(rand($limit));
-  my $end = 1 + int(rand($limit));
+  my $beg = 0 + int(rand($limit));
+  my $end = 0 + int(rand($limit));
   ($beg,$end) = ($end,$beg) if $beg > $end;
 
   # Does range return the same values?

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