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 f400625a52f76dfcafaeda8cb186759eba2b5a11
Author: Dana Jacobsen <d...@acm.org>
Date:   Thu Feb 21 01:08:41 2013 -0800

    Range Moebius function.  Fast Mertens function.
---
 Changes                |  4 ++-
 XS.xs                  | 44 ++++++++++++++++++++++++++++++
 lib/Math/Prime/Util.pm | 74 +++++++++++++++++++++++++++++++++++++++++++++++---
 t/19-moebius.t         | 15 ++++++++--
 util.h                 |  2 ++
 5 files changed, 131 insertions(+), 8 deletions(-)

diff --git a/Changes b/Changes
index 3dfe9c8..dc86582 100644
--- a/Changes
+++ b/Changes
@@ -6,7 +6,9 @@ Revision history for Perl extension Math::Prime::Util.
 
     - Spelling fixes in documentation.
 
-    - Fast return of Euler's totient over a range.
+    - Fast return of Euler's totient and Möbius function over a range.
+
+    - Add mertens function: 500+ times faster than summing moebius($_).
 
     - Speedup of divisor sum.  Also default to sigma if no sub given.
 
diff --git a/XS.xs b/XS.xs
index 36b67ae..0c742c3 100644
--- a/XS.xs
+++ b/XS.xs
@@ -500,3 +500,47 @@ _XS_totient_range(IN UV lo, IN UV hi)
   OUTPUT:
     RETVAL
 
+SV*
+_XS_moebius_range(IN UV lo, IN UV hi)
+  PREINIT:
+    IV* mu;
+    AV* av = newAV();
+    UV i;
+  CODE:
+    /* Return array ref of Moebius function for all n: lo <= n <= hi. */
+    mu = _moebius_range(lo, hi);
+    MPUassert( mu != 0, "_moebius_range returned 0" );
+    for (i = lo; i <= hi; i++)
+      av_push(av,newSViv(mu[i-lo]));
+    Safefree(mu);
+    RETVAL = newRV_noinc( (SV*) av );
+  OUTPUT:
+    RETVAL
+
+IV
+_XS_mertens(IN UV n)
+  PREINIT:
+    IV* mu;
+    IV sum;
+    UV n3, k, startk, endk, limit;
+  CODE:
+    /* Benito and Varona 2008, theorem 3.  Segment. */
+    limit = 1000000;
+    n3 = n/3;
+    sum = 0;
+    startk = 1;
+    endk = limit;
+    prime_precalc( (UV) (sqrt(n)+0.5) );
+    while (startk <= n3) {
+      if (endk > n3) endk = n3;
+      mu = _moebius_range(startk, endk);
+      for (k = startk; k <= endk; k++)
+        if (mu[k-startk] != 0)
+          sum += mu[k-startk] * ((n-k)/(2*k));
+      Safefree(mu);
+      startk = endk+1;
+      endk += limit;
+    }
+    RETVAL = -sum;
+  OUTPUT:
+    RETVAL
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 5300f50..d73e781 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 euler_phi jordan_totient
+                     moebius mertens euler_phi jordan_totient
                      divisor_sum
                      ExponentialIntegral LogarithmicIntegral RiemannZeta 
RiemannR
                    );
@@ -1067,10 +1067,34 @@ sub all_factors {
 # Merten's functions also.
 
 sub moebius {
-  my($n) = @_;
+  my($n, $nend) = @_;
   _validate_positive_integer($n, 1);
-  return 1 if $n == 1;
 
+  # Moebius over a range.
+  if (defined $nend) {
+    _validate_positive_integer($nend);
+    return () if $nend < $n;
+    if ($nend <= $_XS_MAXVAL) {
+      my $mu = _XS_moebius_range($n, $nend);
+      return @$mu;
+    }
+    my @mu = map { 1 } 0 .. $nend;
+    foreach my $j (2 .. $nend) {
+      next unless $mu[$j] == 1;
+      for (my $i = $j; $i <= $nend; $i += $j) {
+        $mu[$i] = ($mu[$i] == 1) ? -$j : -$mu[$i];
+      }
+    }
+    for (my $j = 2; $j*$j <= $nend; $j++) {
+      next unless $mu[$j] == -$j;
+      for (my $i = $j*$j; $i <= $nend; $i += $j*$j) {
+        $mu[$i] = 0;
+      }
+    }
+    return map { ($_>0) - ($_<0) } @mu[ $n .. $nend ];
+  }
+
+  return 1 if $n == 1;
   # Quick check for small replicated factors
   return 0 if ($n >= 25) && (!($n % 4) || !($n % 9) || !($n % 25));
 
@@ -1082,6 +1106,22 @@ sub moebius {
   return (((scalar @factors) % 2) == 0) ? 1 : -1;
 }
 
+sub mertens {
+  my($n) = @_;
+  _validate_positive_integer($n);
+  return (0,1)[$n] if $n <= 1;
+  return _XS_mertens($n) if $n <= $_XS_MAXVAL;
+  # See Benito and Varona 2008, theorem 3
+  my $n3 = int($n/3);
+  my @mu = (0, moebius(1, $n3));
+  my $sum = 0;
+  foreach my $k (1 .. $n3) {
+    next if $mu[$k] == 0;
+    $sum += $mu[$k] * int(($n - $k) / (2*$k));
+  }
+  return -$sum;
+}
+
 
 # Euler Phi, aka Euler Totient.  A000010
 
@@ -1906,6 +1946,8 @@ Version 0.21
 
   # Moebius function used to calculate Mertens
   $sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
+  # Mertens function directly (more efficient for large values)
+  say mertens(10_000_000);
 
   # divisor sum
   $sigma  = divisor_sum( $n );
@@ -2337,6 +2379,30 @@ C<n = 1>, 0 if C<n> is not square free (i.e. C<n> has a 
repeated factor),
 and C<-1^t> if C<n> is a product of C<t> distinct primes.  This is an
 important function in prime number theory.
 
+If called with two arguments, they define a range C<low> to C<high>, and the
+function returns an array with the value of the Möbius function for every n
+from low to high inclusive.  Large values of high will result in a lot of
+memory use.  The algorithm used is Deléglise and Rivat (1996) algorithm 4.1,
+which is a segmented version of Lioen and van de Lune (1994) algorithm 3.2.
+
+
+=head2 mertens
+
+  say "Mertens(10M) = ", mertens(10_000_000);   # = 1037
+
+Returns the Mertens function of the positive non-zero integer input.  This
+function is defined as C<sum(moebius(1..n))>.  This is a much more efficient
+solution for larger inputs.  For example, computing Mertens(100M) takes:
+
+   time    approx mem
+     0.4s      8MB     mertens(100_000_000)
+    74.8s   7000MB     List::Util::sum(moebius(1,100_000_000))
+   325.7s      0MB     $sum += moebius($_) for 1..100_000_000
+
+The algorithm used is a segmented version of the C<n/3> summation from Benito
+and Varona (2008) theorem 3.  A future version may implement an even more
+efficient calculation, such as the algorithm of Deléglise and Rivat (1996).
+
 
 =head2 euler_phi
 
@@ -2350,7 +2416,7 @@ 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,
+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.
 
 
diff --git a/t/19-moebius.t b/t/19-moebius.t
index ed13516..3da2504 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -3,7 +3,7 @@ use strict;
 use warnings;
 
 use Test::More;
-use Math::Prime::Util qw/moebius euler_phi jordan_totient divisor_sum/;
+use Math::Prime::Util qw/moebius mertens euler_phi jordan_totient divisor_sum/;
 
 my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
 my $broken64 = (18446744073709550592 == ~0);
@@ -66,7 +66,7 @@ my %sigmak = (
 
 plan tests => 0 + 1
                 + 1 # Small Moebius
-                + scalar(keys %mertens)
+                + 3*scalar(keys %mertens) + 3
                 + 2 # Small Phi
                 + scalar(keys %totients)
                 + scalar(keys %jordan_totients)
@@ -85,8 +85,17 @@ ok(!eval { moebius(0); }, "moebius(0)");
 while (my($n, $mertens) = each (%mertens)) {
   my $M = 0;
   $M += moebius($_) for (1 .. $n);
-  is( $M, $mertens, "Mertens($n) == $mertens" );
+  is( $M, $mertens, "sum(moebius(k) for k=1..$n) == $mertens" );
+  # Calculate using ranged moebius
+  $M = 0;
+  $M += $_ for moebius(1,$n);
+  is( $M, $mertens, "sum(moebius(1..$n) == $mertens" );
+  # Now with mertens function
+  is( mertens($n), $mertens, "mertens($n) == $mertens" );
 }
+is( mertens(  100000),  -48, "mertens(100000)" );
+is( mertens( 1000000),  212, "mertens(1000000)" );
+is( mertens(10000000), 1037, "mertens(10000000)" );
 
 {
   my @phi = map { euler_phi($_) } (0 .. $#A000010);
diff --git a/util.h b/util.h
index b300a3d..94216f0 100644
--- a/util.h
+++ b/util.h
@@ -15,6 +15,8 @@ extern UV  _XS_prev_prime(UV x);
 extern UV  _XS_prime_count(UV low, UV high);
 extern UV  _XS_nth_prime(UV x);
 
+extern IV*  _moebius_range(UV low, UV high);
+
 extern double _XS_ExponentialIntegral(double x);
 extern double _XS_LogarithmicIntegral(double x);
 extern long double ld_riemann_zeta(long double x);

-- 
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