[libmath-prime-util-perl] 09/17: Range Moebius function. Fast Mertens function.
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:47:14 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 f400625a52f76dfcafaeda8cb186759eba2b5a11
Author: Dana Jacobsen <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list