[libmath-prime-util-perl] 09/23: Add jordan_totient, divisor_sum. Slight speedup for euler_phi and moebius.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:55 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.14
in repository libmath-prime-util-perl.
commit 433250b55a1990d71139077b6a5d3b696a31cb37
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Nov 22 23:35:20 2012 -0800
Add jordan_totient, divisor_sum. Slight speedup for euler_phi and moebius.
---
Changes | 4 +++
TODO | 13 ++++---
lib/Math/Prime/Util.pm | 98 +++++++++++++++++++++++++++++++++++++++++++++++---
t/19-moebius.t | 2 +-
4 files changed, 108 insertions(+), 9 deletions(-)
diff --git a/Changes b/Changes
index e3eae4b..046f700 100644
--- a/Changes
+++ b/Changes
@@ -6,6 +6,10 @@ Revision history for Perl extension Math::Prime::Util.
- XS AKS extended from half-word to full-word.
+ - Add functions:
+ jordan_totient generalization of Euler Totient
+ divisor_sum run coderef for every divisor
+
0.13 19 November 2012
- Fix an issue with prime count, and make prime count available as a
diff --git a/TODO b/TODO
index c35bdee..e317e7f 100644
--- a/TODO
+++ b/TODO
@@ -22,8 +22,6 @@
- After the factoring changes, we need to use Devel::Cover again to ferret
out numbers that pass the early tests.
-- tests for primorial
-
- Test all routines for numbers on word-size boundary, or ranges that cross.
- Test all functions return either native or bigints. Functions that return
@@ -33,8 +31,6 @@
- Add Lehmer in PP (I have a version, just needs some polishing).
-- Move AKS tests into their own test, and add some provable prime tests.
-
- For bignums, RiemannZeta and RiemmannR are slow and give questionable
precision. We should be able to do better. One problem is the accuracy
bug in Math::BigFloat. Perhaps check for Math::MPFR installed and use it?
@@ -42,3 +38,12 @@
- An assembler version of mulmod for i386 would be _really_ helpful for
all the non-x86-64 Intel machines.
+
+- Tests for:
+
+ primorial
+ jordan_totient
+ divisor_sum
+ is_provable_prime
+ RiemannZeta
+ RiemannR
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index a6d1d36..f257e90 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -24,7 +24,9 @@ our @EXPORT_OK = qw(
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime random_nbit_prime random_maurer_prime
primorial pn_primorial
- factor all_factors moebius euler_phi
+ factor all_factors
+ moebius euler_phi jordan_totient
+ divisor_sum
ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
);
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
@@ -737,7 +739,7 @@ sub moebius {
# Quick check for small replicated factors
return 0 if ($n >= 25) && (!($n % 4) || !($n % 9) || !($n % 25));
- my @factors = factor($n);
+ my @factors = ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
my %all_factors;
foreach my $factor (@factors) {
return 0 if $all_factors{$factor}++;
@@ -756,7 +758,8 @@ sub euler_phi {
return 1 if $n <= 1;
my %factor_mult;
- my @factors = grep { !$factor_mult{$_}++ } factor($n);
+ my @factors = grep { !$factor_mult{$_}++ }
+ ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
# Direct from Euler's product formula. Note division will be exact.
#my $totient = $n;
@@ -788,6 +791,54 @@ sub euler_phi {
$totient;
}
+# Jordan's totient -- a generalization of Euler's totient.
+sub jordan_totient {
+ my($k, $n) = @_;
+ _validate_positive_integer($k, 1);
+ return euler_phi($n) if $k == 1;
+
+ return 0 if defined $n && $n <= 0; # Following SAGE's logic here.
+ _validate_positive_integer($n);
+ return 1 if $n <= 1;
+
+ my %factor_mult;
+ my @factors = grep { !$factor_mult{$_}++ }
+ ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
+
+ my $totient = $n - $n + 1;
+
+ if (ref($n) ne 'Math::BigInt') {
+ foreach my $factor (@factors) {
+ my $fmult = int($factor ** $k);
+ $totient *= ($fmult - 1);
+ $totient *= $fmult for (2 .. $factor_mult{$factor});
+ }
+ } else {
+ foreach my $factor (@factors) {
+ my $fmult = $n->copy->bzero->badd("$factor")->bpow($k);
+ $totient->bmul($fmult->copy->bsub(1));
+ $totient->bmul($fmult) for (2 .. $factor_mult{$factor});
+ }
+ }
+ return $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);
+
+ my @afactors = all_factors($n);
+
+ my $sum = $n - $n; # zero as an object of type $n.
+ foreach my $f (1, all_factors($n), $n) {
+ $sum += $sub->($f);
+ }
+ return $sum;
+}
+
# Omega function A001221. Just an example.
sub _omega {
my($n) = @_;
@@ -1479,12 +1530,16 @@ Version 0.13
# Get all factors
@divisors = all_factors( $n );
- # Euler phi (aka the totient) on a large number
+ # Euler phi (Euler's totient) on a large number
use bigint; say euler_phi( 801294088771394680000412 );
+ say jordan_totient(5, 1234); # Jordan's totient
# Moebius function used to calculate Mertens
$sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
+ # divisor sum
+ $sigma = divisor_sum( $n, sub { $_[0] } );
+
# The primorial n# (product of all primes <= n)
say "15# (2*3*5*7*11*13) is ", primorial(15);
# The primorial p(n)# (product of first n primes)
@@ -1912,6 +1967,41 @@ 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>.
+=head2 jordan_totient
+
+ say "Jordan's totient J_$k($n) is ", jordan_totient($k, $n);
+
+Returns Jordan's totient function for a given integer value. Jordan's totient
+is a generalization of Euler's totient, where
+ C<jordan_totient(1,$n) == euler_totient($n)>
+This counts the number of k-tuples less than or equal to n that form a coprime
+tuple with n. As with C<euler_phi>, 0 is returned for all C<n E<lt> 1>.
+This function can be used to generate some other useful functions, such as
+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] } );
+
+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:
+
+ divisor_sum( $n, sub { my $d=shift; $d**5 * moebius($n/$d); } );
+
+calculates the 5th Jordan totient (OEIS 59378). In this example we have a
+specific function C<jordan_totient> that can compute this more efficiently.
+
+ 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.
+
+
=head2 primorial
$prim = primorial(11); # 11# = 2*3*5*7*11 = 2310
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 0500d35..55b1caa 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -16,7 +16,7 @@ my %mertens = (
100 => 1,
1000 => 2,
10000 => -23,
- 100000 => -48,
+# 100000 => -48,
# 1000000 => 212,
# 10000000 => 1037,
);
--
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