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 <d...@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

_______________________________________________
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