[libmath-prime-util-perl] 29/72: Divisor sum with integer second arg
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:38 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list