[libmath-prime-util-perl] 19/72: Add speedup for divisor count

[Partha P. Mukherjee]

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

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

commit c2a0056165a28e490109adf3e4d907d712990c2d
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Sep 5 13:50:30 2013 -0700

    Add speedup for divisor count
---
 Changes                |  3 +++
 lib/Math/Prime/Util.pm | 21 +++++++++++++++++++++
 t/19-moebius.t         | 10 +++++++++-
 3 files changed, 33 insertions(+), 1 deletion(-)

diff --git a/Changes b/Changes
index 499227f..695ed5b 100644
--- a/Changes
+++ b/Changes
@@ -17,6 +17,9 @@ Revision history for Perl module Math::Prime::Util
       should stay under ~150MB even with giant sizes.
       Thanks to Kim Walisch for all discussions about this.
 
+    - divisor_sum can take a '1' as the second argument as an alias for
+      sub {1}, but works much faster.
+
 0.31  2013-08-07
 
     - Change proof certificate documentation to reflect the new text format.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 7399339..a2d7743 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1437,6 +1437,22 @@ sub divisor_sum {
     return $product;
   }
 
+  if (ref($sub) ne 'CODE' && int($sub) == 1) {
+    return 1 if $n == 1;
+    my ($product, $exponent) = (1, 1);
+    my @factors = ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
+    while (@factors) {
+      if (@factors > 1 && $factors[0] == $factors[1]) {
+        $exponent++;
+      } else {
+        $product *= ($exponent+1);
+        $exponent = 1;
+      }
+      shift @factors;
+    }
+    return $product;
+  }
+
   croak "Second argument must be a code ref" unless ref($sub) eq 'CODE';
   my $sum = $sub->(1);
   return $sum if $n == 1;
@@ -2419,6 +2435,7 @@ Version 0.31
 
   # divisor sum
   $sigma  = divisor_sum( $n );
+  $sigma0 = divisor_sum( $n, 1 );
   $sigma2 = divisor_sum( $n, sub { $_[0]*$_[0] } );
 
   # primorial n#, primorial p(n)#, and lcm
@@ -3394,6 +3411,10 @@ This function takes a positive integer as input and returns the sum of all
 the divisors of the input, including 1 and itself.  This is known as the
 sigma function (see Hardy and Wright section 16.7, or OEIS A000203).
 
+If a second argument is given that is numerically equal to 1, then we do
+a summation of the number of divisors.  This is identical to the general
+form below with C<sub { 1 }>, but runs faster.
+
 The more general form takes a code reference as a second parameter, which
 is applied to each divisor before the summation.  This allows computation
 of numerous functions such as OEIS A000005 [d(n), sigma_0(n), tau(n)]:
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 057f8ff..2e3796e 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -169,7 +169,7 @@ plan tests => 0 + 1
                 + 2  # Dedekind psi calculated two ways
                 + 1  # Calculate J5 two different ways
                 + 2 * $use64 # Jordan totient example
-                + 1 + scalar(keys %sigmak)
+                + 1 + scalar(keys %sigmak) + 2
                 + scalar(keys %mangoldt)
                 + scalar(keys %chebyshev1)
                 + scalar(keys %chebyshev2);
@@ -263,6 +263,14 @@ while (my($k, $sigmaref) = each (%sigmak)) {
   my @slist = map { divisor_sum($_) } 1 .. scalar @{$sigmak{1}};
   is_deeply(\@slist, $sigmak{1}, "divisor_sum(n)");
 }
+# tau two ways
+{
+  my $len = scalar @{$sigmak{0}};
+  my @slist1 = map { divisor_sum($_, sub {1}) } 1 .. $len;
+  my @slist2 = map { divisor_sum($_, 1      ) } 1 .. $len;
+  is_deeply( \@slist1, $sigmak{0}, "tau as divisor_sum(n, sub {1})" );
+  is_deeply( \@slist2, $sigmak{0}, "tau as divisor_sum(n, 1)" );
+}
 
 ###### Exponential of von Mangoldt
 while (my($n, $em) = each (%mangoldt)) {

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