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 9cd55d37ddabe5aa40841bfe10be1f2e9382dc3e
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Sep 13 19:25:21 2013 -0700

    Divisor sum with integer second arg
---
 Changes                |  5 +--
 XS.xs                  |  5 +--
 factor.c               | 38 +++++++++++++++++------
 factor.h               |  2 +-
 lib/Math/Prime/Util.pm | 82 ++++++++++++++++++++++++++------------------------
 t/19-moebius.t         | 11 +++++--
 6 files changed, 86 insertions(+), 57 deletions(-)

diff --git a/Changes b/Changes
index cc86e62..d41e4f0 100644
--- a/Changes
+++ b/Changes
@@ -18,8 +18,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.  (experimental, may change)
+    - divisor_sum can take an integer 'k' in the second argument to compute
+      sigma_k.  This is much faster than using subs, especially when the
+      result can be computed in XS using native precision.
 
     - Incorporate Montgomery reduction for primality testing, thanks to
       Wojciech Izykowski.  This is a 1.3 to 1.5x speedup for is_prob_prime,
diff --git a/XS.xs b/XS.xs
index 84a7ba3..99c0e64 100644
--- a/XS.xs
+++ b/XS.xs
@@ -188,7 +188,6 @@ _XS_nth_prime(IN UV n)
     _XS_meissel_pi = 4
     _XS_lehmer_pi = 5
     _XS_LMO_pi = 6
-    _XS_divisor_sum = 7
   PREINIT:
     UV ret;
   CODE:
@@ -200,13 +199,15 @@ _XS_nth_prime(IN UV n)
       case 4: ret = _XS_meissel_pi(n); break;
       case 5: ret = _XS_lehmer_pi(n); break;
       case 6: ret = _XS_LMO_pi(n); break;
-      case 7: ret = _XS_divisor_sum(n); break;
       default: croak("_XS_nth_prime: Unknown function alias"); break;
     }
     RETVAL = ret;
   OUTPUT:
     RETVAL
 
+UV
+_XS_divisor_sum(IN UV n, IN UV k)
+
 
 UV
 _get_prime_cache_size()
diff --git a/factor.c b/factor.c
index e43f335..94e919c 100644
--- a/factor.c
+++ b/factor.c
@@ -291,7 +291,7 @@ static int is_perfect_square(UV n, UV* sqrtn)
 }
 
 
-UV _XS_divisor_sum(UV n)
+UV _XS_divisor_sum(UV n, UV k)
 {
   UV factors[MPU_MAX_FACTORS+1];
   int nfac, i;
@@ -300,14 +300,34 @@ UV _XS_divisor_sum(UV n)
   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;
+    UV f = factors[i];
+    UV e = 1;
+    while (i+1 < nfac && f == factors[i+1])
+      { e++; i++; }
+    if (k == 0) {
+      product *= (e+1);
+    } else if (k == 1 && e == 1) {
+      product *= f+1;
+    } else if (k == 1) {
+      UV fmult = f * f;
+      while (e-- >= 2)
+        fmult *= f;
+      product *= (fmult - 1) / (f - 1);
+    } else {   /* Overflow is a concern */
+      UV j, pk = f * f;
+      for (j = 2; j < k; j++)  pk *= f;
+      if (e == 1) {
+        product *= pk+1;
+      } else {
+        /* Less overflow than (pk^(e+1)-1) / (pk-1) */
+        UV fmult = 1 + pk + pk*pk;
+        UV pke = pk * pk;
+        for (j = 2; j < e; j++) {
+          pke *= pk;
+          fmult += pke;
+        }
+        product *= fmult;
+      }
     }
   }
   return product;
diff --git a/factor.h b/factor.h
index 2979b32..06ea83d 100644
--- a/factor.h
+++ b/factor.h
@@ -18,6 +18,6 @@ extern int pplus1_factor(UV n, UV *factors, UV B);
 extern int squfof_factor(UV n, UV *factors, UV rounds);
 extern int racing_squfof_factor(UV n, UV *factors, UV rounds);
 
-extern UV _XS_divisor_sum(UV n);
+extern UV _XS_divisor_sum(UV n, UV k);
 
 #endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 556787f..533ba76 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1413,52 +1413,53 @@ sub jordan_totient {
 
 # Mathematica and Pari both have functions like this.
 sub divisor_sum {
-  my($n, $sub) = @_;
+  my($n, $k) = @_;
   return (0,1)[$n] if defined $n && $n <= 1;
   _validate_num($n) || _validate_positive_integer($n);
 
-  if (!defined $sub) {
-    return _XS_divisor_sum($n) if $n <= $_XS_MAXVAL;
-    my $bone = (ref($n) eq 'Math::BigInt') ? $n->copy->bone : 1;
-    my $product = $bone;
-    my @factors = factor($n);
-    while (@factors) {
-      if (@factors > 1 && $factors[0] == $factors[1]) {
-        my $fmult = $bone * $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;
+  $k = 1 unless defined $k;
+
+  if (ref($k) eq 'CODE') {
+    my $sum = $k->(1);
+    return $sum if $n == 1;
+    foreach my $f (all_factors($n), $n ) {
+      $sum += $k->($f);
     }
-    return $product;
+    return $sum;
   }
 
-  # TODO: Alternately the argument could be the k for sigma, so this
-  # should be 0, the above should be 1, etc.
-  if (ref($sub) ne 'CODE' && int($sub) == 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 or number"
+    unless _validate_num($k) || _validate_positive_integer($k);
+
+  if ($n <= $_XS_MAXVAL) {
+    return _XS_divisor_sum($n, $k) if $k <= 1;
+    # 2 overflows 64-bit at 3,000,000,000 ish
+    # 3 overflows 64-bit at 2487240
+    # 4 overflows 64-bit at 64260
+    # If we think we can handle overflow, we could do a few more
   }
 
-  croak "Second argument must be a code ref" unless ref($sub) eq 'CODE';
-  my $sum = $sub->(1);
-  return $sum if $n == 1;
-  foreach my $f (all_factors($n), $n ) {
-    $sum += $sub->($f);
+  my $bone = (ref($_[0]) eq 'Math::BigInt') ? $_[0]->copy->bone : 1;
+  my $product = $bone;
+  my @factors = ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
+  while (@factors) {
+    my $f = shift @factors;
+    my $e = 1;
+    while (@factors && $factors[0] == $f) {
+      $e++;
+      shift @factors;
+    }
+    if ($k == 0)         { $product *= ($e+1); }
+    elsif ($k == 1) {
+      if ($e == 1)       { $product *= $f+1; }
+      else               { $product *= ( ($bone*$f)**($e+1) - 1) / ($f - 1); }
+    } else {
+      # TODO: do the long way around to get less overflow, see factor.c
+      my $pk = ($bone * $f) ** $k;
+      $product *= ($e == 1)  ?  $pk+1  :  ($pk ** ($e+1) - 1) / ($pk - 1);
+    }
   }
-  return $sum;
+  return $product;
 }
 
                                    # Need proto for the block
@@ -3453,9 +3454,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.
+If a second argument is given that is numeric, it is used as the sigma
+parameter.  The result will then be the sum of each divisor raised to
+this power.  Hence C<0> will count the divisors, C<1> will sum them,
+C<2> will sum the squares, and so on.
 
 The more general form takes a code reference as a second parameter, which
 is applied to each divisor before the summation.  This allows computation
diff --git a/t/19-moebius.t b/t/19-moebius.t
index b44b888..dbc75ef 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -194,7 +194,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) + 2
+                + 1 + 2*scalar(keys %sigmak) + 2
                 + scalar(keys %mangoldt)
                 + scalar(keys %chebyshev1)
                 + scalar(keys %chebyshev2);
@@ -282,6 +282,11 @@ while (my($k, $sigmaref) = each (%sigmak)) {
     push @slist, divisor_sum( $n, sub { int($_[0] ** $k) } );
   }
   is_deeply( \@slist, $sigmaref, "Sum of divisors to the ${k}th power: 
Sigma_$k" );
+  @slist = ();
+  foreach my $n (1 .. scalar @$sigmaref) {
+    push @slist, divisor_sum( $n, $k );
+  }
+  is_deeply( \@slist, $sigmaref, "Sigma_$k using integer instead of sub" );
 }
 # k=1 standard sum -- much faster
 {
@@ -292,9 +297,9 @@ while (my($k, $sigmaref) = each (%sigmak)) {
 {
   my $len = scalar @{$sigmak{0}};
   my @slist1 = map { divisor_sum($_, sub {1}) } 1 .. $len;
-  my @slist2 = map { divisor_sum($_, 1      ) } 1 .. $len;
+  my @slist2 = map { divisor_sum($_, 0      ) } 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)" );
+  is_deeply( \@slist2, $sigmak{0}, "tau as divisor_sum(n, 0)" );
 }
 
 ###### Exponential of von Mangoldt

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

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to