[libmath-prime-util-perl] 06/17: Add Euler totient range, speedup divisor_sum
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:47:13 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.21
in repository libmath-prime-util-perl.
commit ccc7965299864634b3928539a659473807f41ae8
Author: Dana Jacobsen <dana at acm.org>
Date: Tue Feb 19 04:00:22 2013 -0800
Add Euler totient range, speedup divisor_sum
---
Changes | 4 ++
XS.xs | 39 +++++++++++++
lib/Math/Prime/Util.pm | 150 ++++++++++++++++++++++++++++++++++++-------------
t/19-moebius.t | 6 +-
4 files changed, 158 insertions(+), 41 deletions(-)
diff --git a/Changes b/Changes
index c788a06..46ec96c 100644
--- a/Changes
+++ b/Changes
@@ -6,6 +6,10 @@ Revision history for Perl extension Math::Prime::Util.
- Spelling fixes in documentation.
+ - Fast return of Euler's totient over a range.
+
+ - Speedup of divisor sum. Also default to sigma if no sub given.
+
0.20 3 February 2012
- Speedup for PP AKS, and turn off test on 32-bit machines.
diff --git a/XS.xs b/XS.xs
index 15f51ff..36b67ae 100644
--- a/XS.xs
+++ b/XS.xs
@@ -461,3 +461,42 @@ _XS_RiemannZeta(double x)
double
_XS_RiemannR(double x)
+
+
+SV*
+_XS_totient_range(IN UV lo, IN UV hi)
+ PREINIT:
+ UV* totients;
+ AV* av = newAV();
+ UV i, j;
+ CODE:
+ /* Calculate Euler's totient for all n: lo <= n <= hi. */
+ /* Return as array ref */
+ New(0, totients, hi+1, UV);
+ if (totients == 0)
+ croak("Could not get memory for %"UVuf" totients\n", hi);
+ for (i = 0; i <= hi; i++)
+ totients[i] = i;
+ if (lo <= 0 && hi >= 0) av_push(av,newSVuv(totients[0]));
+ if (lo <= 1 && hi >= 1) av_push(av,newSVuv(totients[1]));
+ if (lo <= 2 && hi >= 2) {
+ totients[2] = 1;
+ av_push(av,newSVuv(totients[2]));
+ }
+ for (j = 2; j <= hi/2; j++)
+ totients[2*j] /= 2;
+ for (i = 3; i <= hi; i++) {
+ if (totients[i] == i) {
+ totients[i] = i-1;
+ for (j = 2; j <= hi/i; j++) {
+ totients[i*j] = (totients[i*j]*(i-1))/i;
+ }
+ }
+ if (i >= lo)
+ av_push(av,newSVuv(totients[i]));
+ }
+ Safefree(totients);
+ RETVAL = newRV_noinc( (SV*) av );
+ OUTPUT:
+ RETVAL
+
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index f40d403..ea3f0b1 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1038,26 +1038,28 @@ sub pn_primorial {
sub all_factors {
my $n = shift;
+ my $use_bigint = defined $bigint::VERSION
+ || !($n < $_Config{'maxparam'} || int($n) eq $_Config{'maxparam'});
my @factors = factor($n);
my %all_factors;
- foreach my $f1 (@factors) {
- next if $f1 >= $n;
- # We're adding to %all_factors in the loop, so grab the keys now.
- my @all = keys %all_factors;;
- if (!defined $bigint::VERSION) {
- foreach my $f2 (@all) {
- $all_factors{$f1*$f2} = 1 if ($f1*$f2) < $n;
- }
- } else {
- # Many of the factors will be numified after coming back, so we need
- # to make sure we're using bigints when we calculate the product.
- foreach my $f2 (@all) {
- my $product = Math::BigInt->new("$f1") * Math::BigInt->new("$f2");
- $product = int($product->bstr) if $product <= ~0;
- $all_factors{$product} = 1 if $product < $n;
- }
+ if ($use_bigint) {
+ foreach my $f1 (@factors) {
+ next if $f1 >= $n;
+ my $big_f1 = Math::BigInt->new("$f1");
+ my @to_add = map { ($_ <= ~0) ? int($_->bstr) : $_ }
+ grep { $_ < $n }
+ map { $big_f1 * $_ }
+ keys %all_factors;
+ undef @all_factors{ $f1, @to_add };
+ }
+ } else {
+ foreach my $f1 (@factors) {
+ next if $f1 >= $n;
+ my @to_add = grep { $_ < $n }
+ map { int($f1 * $_) }
+ keys %all_factors;
+ undef @all_factors{ $f1, @to_add };
}
- $all_factors{$f1} = 1;
}
@factors = sort {$a<=>$b} keys %all_factors;
return @factors;
@@ -1090,12 +1092,35 @@ sub moebius {
# Euler Phi, aka Euler Totient. A000010
sub euler_phi {
- my($n) = @_;
+ my($n, $nend) = @_;
# SAGE defines this to be 0 for all n <= 0. Others choose differently.
- return 0 if defined $n && $n <= 0; # Following SAGE's logic here.
+ # I am following SAGE's decision for n <= 0.
+ return 0 if defined $n && $n < 0;
_validate_positive_integer($n);
- return 1 if $n <= 1;
+ undef $nend if defined $nend && $nend == $n;
+
+ # Totient over a range. Could be improved, as this can use huge memory.
+ if (defined $nend) {
+ _validate_positive_integer($nend);
+ return () if $nend < $n;
+ # Use XS code if at all possible. Better memory use.
+ if ($nend <= $_XS_MAXVAL) {
+ my $totients = _XS_totient_range($n, $nend);
+ return @$totients;
+ }
+ my @totients = map { $_ } 0 .. $nend;
+ foreach my $i (2 .. $nend) {
+ next unless $totients[$i] == $i;
+ $totients[$i] = $i-1;
+ foreach my $j (2 .. int($nend / $i)) {
+ $totients[$i*$j] = ($totients[$i*$j]*($i-1))/$i;
+ }
+ }
+ splice(@totients, 0, $n) if $n > 0;
+ return @totients;
+ }
+ return (0,1)[$n] if $n <= 1;
my %factor_mult;
my @factors = grep { !$factor_mult{$_}++ }
($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
@@ -1153,8 +1178,9 @@ sub jordan_totient {
$totient *= $fmult for (2 .. $factor_mult{$factor});
}
} else {
+ my $zero = $n->copy->bzero;
foreach my $factor (@factors) {
- my $fmult = $n->copy->bzero->badd("$factor")->bpow($k);
+ my $fmult = $zero->copy->badd("$factor")->bpow($k);
$totient->bmul($fmult->copy->bsub(1));
$totient->bmul($fmult) for (2 .. $factor_mult{$factor});
}
@@ -1165,15 +1191,23 @@ sub jordan_totient {
# Mathematica and Pari both have functions like this.
sub divisor_sum {
my($n, $sub) = @_;
- croak "Second argument must be a code ref" unless ref($sub) eq 'CODE';
return 0 if defined $n && $n < 1;
_validate_positive_integer($n);
- return ($n-$n+$sub->(1)) if $n == 1;
- my @afactors = all_factors($n);
+ if (!defined $sub) {
+ return $n if $n == 1;
+ my $sum = 1 + $n;
+ foreach my $f ( all_factors($n) ) {
+ $sum += $f;
+ }
+ return $sum;
+ }
- my $sum = $n - $n; # zero as an object of type $n.
- foreach my $f (1, all_factors($n), $n) {
+ croak "Second argument must be a code ref" unless ref($sub) eq 'CODE';
+ my $sum = $n - $n;
+ return ($sum+$sub->(1)) if $n == 1;
+
+ foreach my $f (1, all_factors($n), $n ) {
$sum += $sub->($f);
}
return $sum;
@@ -1880,7 +1914,8 @@ Version 0.21
$sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
# divisor sum
- $sigma = divisor_sum( $n, sub { $_[0] } );
+ $sigma = divisor_sum( $n );
+ $sigma2 = divisor_sum( $n, sub { $_[0]*$_[0] } );
# The primorial n# (product of all primes <= n)
say "15# (2*3*5*7*11*13) is ", primorial(15);
@@ -2320,6 +2355,10 @@ C<n>. Given the definition used, C<euler_phi> will return 0 for all
C<n E<lt> 1>. This follows the logic used by SAGE. Mathematic/WolframAlpha
also returns 0 for input 0, but returns C<euler_phi(-n)> for C<n E<lt> 0>.
+If called with two arguments, they define a range C<low> to C<high>, and the
+function returns an array with the totient of every n from low to high,
+inclusive. Large values of high will result in a lot of memory use.
+
=head2 jordan_totient
@@ -2336,24 +2375,38 @@ the Dedikind psi function, where C<psi(n) = J(2,n) / J(1,n)>.
=head2 divisor_sum
- say "Sum of divisors of $n:", divisor_sum( $n, sub { $_[0] } );
+ say "Sum of divisors of $n:", divisor_sum( $n );
-This function takes a positive integer as input, along with a code reference.
-For each positive divisor of the input, including 1 and itself, the coderef
-is called with the divisor as the only argument, and the return values summed.
-There are a number of utilities this can be used for, though it may not be the
-most efficient way to calculate them. Example:
+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).
- divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); } );
+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)]:
+
+ divisor_sum( $n, sub { 1 } );
-calculates the 5th Jordan totient (OEIS 59378). In this example we have a
-specific function C<jordan_totient> that can compute this more efficiently.
+OEIS A001157 [sigma_2(n)]:
+
+ divisor_sum( $n, sub { $_[0]*$_[0] } )
+
+the general sigma_k (OEIS A00005, A000203, A001157, A001158, etc.):
divisor_sum( $n, sub { $_[0] ** $k } );
-calculates sigma_k (OEIS A000005, A000203, A001157, A001158 for k=0..3). The
-simple sigma shown as the first example can be used to find aliquot sums,
-abundant numbers, perfect numbers, and more.
+the 5th Jordan totient (OEIS A059378):
+
+ divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); } );
+
+though in the last case we have a function L</"jordan_totient"> to compute
+it more efficiently.
+
+This function is useful for calculating things like aliquot sums, abundant
+numbers, perfect numbers, etc.
+
+The summation is done as a bigint if the input was a bigint object. You may
+need to ensure the result of the subroutine does not overflow a native int.
=head2 primorial
@@ -2941,7 +2994,7 @@ Project Euler, problem 10 (summation of primes):
Project Euler, problem 21 (Amicable numbers):
use Math::Prime::Util qw/divisor_sum/;
- sub dsum { my $n = shift; divisor_sum($n, sub {$_[0]}) - $n; }
+ sub dsum { my $n = shift; divisor_sum($n) - $n; }
my $sum = 0;
foreach my $a (1..10000) {
my $b = dsum($a);
@@ -2979,6 +3032,23 @@ Here's the right way to do PE problem 69 (under 0.03s):
$n++ while pn_primorial($n+1) < 1000000;
say pn_primorial($n);'
+Project Euler, problem 187, stupid brute force solution:
+
+ use Math::Prime::Util qw/factor -nobigint/;
+ my $nsemis = 0;
+ do { my @f = factor($_); $nsemis++ if scalar @f == 2; }
+ for 1 .. int(10**8)-1;
+ say $nsemis;
+
+Here's the best way for PE187. Under 30 milliseconds from the command line:
+
+ use Math::Prime::Util qw/primes prime_count -nobigint/;
+ use List::Util qw/sum/;
+ my $limit = shift || int(10**8);
+ my @primes = @{primes(int(sqrt($limit)))};
+ say sum( map { prime_count(int(($limit-1)/$primes[$_-1])) - $_ + 1 }
+ 1 .. scalar @primes );
+
=head1 LIMITATIONS
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 56b0eb0..ed13516 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -67,7 +67,7 @@ my %sigmak = (
plan tests => 0 + 1
+ 1 # Small Moebius
+ scalar(keys %mertens)
- + 1 # Small Phi
+ + 2 # Small Phi
+ scalar(keys %totients)
+ scalar(keys %jordan_totients)
+ 2 # Dedekind psi calculated two ways
@@ -92,6 +92,10 @@ while (my($n, $mertens) = each (%mertens)) {
my @phi = map { euler_phi($_) } (0 .. $#A000010);
is_deeply( \@phi, \@A000010, "euler_phi 0 .. $#A000010" );
}
+{
+ my @phi = euler_phi(0, $#A000010);
+ is_deeply( \@phi, \@A000010, "euler_phi with range: 0, $#A000010" );
+}
while (my($n, $phi) = each (%totients)) {
is( euler_phi($n), $phi, "euler_phi($n) == $phi" );
}
--
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