[libmath-prime-util-perl] 13/17: Add Mangoldt function

[Partha P. Mukherjee]

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 5c7004570b2210d32d05dff2892d5e59281ab31d
Author: Dana Jacobsen <dana at acm.org>
Date:   Fri Feb 22 09:20:41 2013 -0800

    Add Mangoldt function
---
 Changes                |  4 +++-
 lib/Math/Prime/Util.pm | 35 ++++++++++++++++++++++++++++++++++-
 util.c                 |  1 +
 3 files changed, 38 insertions(+), 2 deletions(-)

diff --git a/Changes b/Changes
index 0623989..714a7af 100644
--- a/Changes
+++ b/Changes
@@ -12,7 +12,9 @@ Revision history for Perl extension Math::Prime::Util.
 
     - Add mertens function: 1000+ times faster than summing moebius($_).
 
-    - Speedup of divisor sum.  Also default to sigma if no sub given.
+    - Add exp_mangoldt function: exponential of von Mangoldt's function.
+
+    - divisor sum is 2x faster.  Also default to sigma if no sub given.
 
 0.20  3 February 2012
 
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 3b4f4ed..3d111e2 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -27,7 +27,7 @@ our @EXPORT_OK = qw(
                      random_strong_prime random_maurer_prime
                      primorial pn_primorial
                      factor all_factors
-                     moebius mertens euler_phi jordan_totient
+                     moebius mertens euler_phi jordan_totient exp_mangoldt
                      divisor_sum
                      ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
                    );
@@ -1271,6 +1271,24 @@ sub _omega {
   return scalar @factors;
 }
 
+# Exponential of Mangoldt function (A014963).
+# Return p if n = p^m [p prime, m >= 1], 1 otherwise.
+sub exp_mangoldt {
+  my($n) = @_;
+  return 1 if defined $n && $n <= 1;
+  _validate_positive_integer($n);
+
+  #my $is_prime = ($n<=$_XS_MAXVAL) ? _XS_is_prob_prime($n) : is_prob_prime($n);
+  #return $n if $is_prime;
+
+  my %factor_mult;
+  my @factors = grep { !$factor_mult{$_}++ }
+                ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
+
+  return 1 unless scalar @factors == 1;
+  return $factors[0];
+}
+
 
 #############################################################################
 # Front ends to functions.
@@ -1963,6 +1981,9 @@ Version 0.21
   # Mertens function directly (more efficient for large values)
   say mertens(10_000_000);
 
+  # Exponential of Mangoldt function
+  say "lamba(49) = ", log(exp_mangoldt(49));
+
   # divisor sum
   $sigma  = divisor_sum( $n );
   $sigma2 = divisor_sum( $n, sub { $_[0]*$_[0] } );
@@ -2453,6 +2474,18 @@ 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 exp_mangoldt
+
+  say "exp(lambda($_)) = ", exp_mangoldt($_) for 1 .. 100;
+
+The Mangoldt function Λ(n) (also known as von Mangoldt's function) is equal
+to log p if n is prime or a power of a prime, and 0 otherwise.  We return
+the exponential so all results are integers.  Hence the return value
+for C<exp_mangoldt> is:
+   C<p> if C<n = p^m> for some prime C<p> and integer C<m E<gt>= 1>
+   1 otherwise.
+
+
 =head2 divisor_sum
 
   say "Sum of divisors of $n:", divisor_sum( $n );
diff --git a/util.c b/util.c
index c09d49f..402106b 100644
--- a/util.c
+++ b/util.c
@@ -642,6 +642,7 @@ IV* _moebius_range(UV lo, UV hi)
    */
   range = hi-lo+1;
   sqrtn = (UV) (sqrt(hi) + 0.5);
+  /* sqrtn = (UV) sqrt(hi);  if (sqrtn*sqrtn < hi) sqrtn++; */  /* TODO */
 
   New(0, mu, range, IV);
   if (mu == 0)

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