[libmath-prime-util-perl] 11/20: 10x speedup for divisor_sum

[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 8aea368b9e763da2630463f2fa5812ea21485d76
Author: Dana Jacobsen <dana at acm.org>
Date:   Sat Mar 2 01:16:33 2013 -0800

    10x speedup for divisor_sum
---
 Changes                |  2 ++
 XS.xs                  |  3 +++
 factor.c               | 24 ++++++++++++++++++++++++
 factor.h               |  2 ++
 lib/Math/Prime/Util.pm | 37 +++++++++++++++++++++++++++++++------
 5 files changed, 62 insertions(+), 6 deletions(-)

diff --git a/Changes b/Changes
index e52cf01..4d3c052 100644
--- a/Changes
+++ b/Changes
@@ -20,6 +20,8 @@ Revision history for Perl extension Math::Prime::Util.
 
     - Add the first and second Chebyshev functions (theta and psi).
 
+    - Divisor sum with no sub is ~10x faster.
+
 0.22 26 February 2013
 
     - Move main factor loop out of xs and into factor.c.
diff --git a/XS.xs b/XS.xs
index 83fc0ce..490b880 100644
--- a/XS.xs
+++ b/XS.xs
@@ -469,3 +469,6 @@ _XS_chebyshev_theta(IN UV n)
 
 double
 _XS_chebyshev_psi(IN UV n)
+
+UV
+_XS_divisor_sum(IN UV n)
diff --git a/factor.c b/factor.c
index b86b02b..8bc1357 100644
--- a/factor.c
+++ b/factor.c
@@ -424,6 +424,30 @@ int _XS_is_prob_prime(UV n)
 }
 
 
+UV _XS_divisor_sum(UV n)
+{
+  UV factors[MPU_MAX_FACTORS+1];
+  int nfac, i;
+  UV product = 1;
+
+  if (n <= 1) return n;
+  nfac = factor(n, factors);
+  for (i = 0; i < nfac; i++) {
+    if (i+1 < nfac && factors[i] == factors[i+1]) {
+      UV fmult = factors[i]*factors[i];
+      do {
+        fmult *= factors[i++];
+      } while (i+1 < nfac && factors[i] == factors[i+1]);
+      product *= (fmult-1) / (factors[i]-1);
+    } else {
+      product *= factors[i]+1;
+    }
+  }
+  return product;
+}
+
+
+
 
 /* Knuth volume 2, algorithm C.
  * Very fast for small numbers, grows rapidly.
diff --git a/factor.h b/factor.h
index b5cf735..a13b52e 100644
--- a/factor.h
+++ b/factor.h
@@ -25,6 +25,8 @@ extern int pminus1_factor(UV n, UV *factors, UV B1, UV B2);
 extern int _XS_miller_rabin(UV n, const UV *bases, int nbases);
 extern int _XS_is_prob_prime(UV n);
 
+extern UV _XS_divisor_sum(UV n);
+
 static UV gcd_ui(UV x, UV y) {
   UV t;
   if (y < x) { t = x; x = y; y = t; }
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index cc7129e..e491903 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1289,16 +1289,26 @@ sub jordan_totient {
 # Mathematica and Pari both have functions like this.
 sub divisor_sum {
   my($n, $sub) = @_;
-  return 0 if defined $n && $n < 1;
+  # I really need to get cracking on an XS validator.
+  #return _XS_divisor_sum($n) if !defined $sub && defined $n && $n <= $_XS_MAXVAL && $_Config{'nobigint'};
+  return (0,1)[$n] if defined $n && $n <= 1;
   _validate_positive_integer($n);
 
   if (!defined $sub) {
-    return $n if $n == 1;
-    my $sum = 1 + $n;
-    foreach my $f ( all_factors($n) ) {
-      $sum += $f;
+    return _XS_divisor_sum($n) if $n <= $_XS_MAXVAL;
+    my $product = 1;
+    my @factors = factor($n);
+    while (@factors) {
+      if (@factors > 1 && $factors[0] == $factors[1]) {
+        my $fmult = $factors[0] * $factors[0];
+        $fmult *= shift @factors while @factors > 1 && $factors[0] == $factors[1];
+        $product *= ($fmult -1) / ($factors[0] - 1);
+      } else {
+        $product *= $factors[0]+1;
+      }
+      shift @factors;
     }
-    return $sum;
+    return $product;
   }
 
   croak "Second argument must be a code ref" unless ref($sub) eq 'CODE';
@@ -3191,6 +3201,21 @@ Print some primes above 64-bit range:
     # Similar code using Pari:
     # perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }'
 
+Examining the η3(x) function of Planat and Solé (2011):
+
+  sub nu3 {
+    my $n = shift;
+    my $phix = chebyshev_psi($n);
+    my $nu3 = 0;
+    foreach my $nu (1..3) {
+      $nu3 += (moebius($nu)/$nu)*LogarithmicIntegral($phix**(1/$nu));
+    }
+    return $nu3;
+  }
+  say prime_count(1000000);
+  say prime_count_approx(1000000);
+  say nu3(1000000);
+
 Project Euler, problem 3 (Largest prime factor):
 
   use Math::Prime::Util qw/factor/;

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