[libmath-prime-util-perl] 01/20: Optimizations for von Mangoldt function
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:47:30 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.23
in repository libmath-prime-util-perl.
commit 3bf8b1e3f094784a95a2601dded06d40e521d783
Author: Dana Jacobsen <dana at acm.org>
Date: Tue Feb 26 13:43:40 2013 -0800
Optimizations for von Mangoldt function
---
Changes | 9 +++++++--
XS.xs | 27 +++++++++++++++++++++++++++
lib/Math/Prime/Util.pm | 46 +++++++++++++++++++++++++---------------------
t/19-moebius.t | 34 ++++++++++++++++++++++++++++++++--
4 files changed, 91 insertions(+), 25 deletions(-)
diff --git a/Changes b/Changes
index d8e8a05..46435b8 100644
--- a/Changes
+++ b/Changes
@@ -1,5 +1,9 @@
Revision history for Perl extension Math::Prime::Util.
+0.23 xx February 2012
+
+ - exp_mangoldt in XS, and speed up the PP version.
+
0.22 26 February 2012
- Move main factor loop out of xs and into factor.c.
@@ -41,8 +45,9 @@ Revision history for Perl extension Math::Prime::Util.
- Speedup for PP AKS, and turn off test on 32-bit machines.
- Replaced fast sqrt detection in PP.pm with a slightly slower version.
- The bloom filter doesn't work right in 32-bit Perl. These two changes
- should speed up testing on some machines by a huge amount.
+ The bloom filter doesn't work right in 32-bit Perl. Having a non-working
+ detector led to really bad performance. Hence this and the AKS change
+ should speed up testing on some 32-bit machines by a huge amount.
- Fix is_perfect_power in XS AKS.
diff --git a/XS.xs b/XS.xs
index fca7d50..95e4733 100644
--- a/XS.xs
+++ b/XS.xs
@@ -3,6 +3,10 @@
#include "perl.h"
#include "XSUB.h"
/* We're not using anything for which we need ppport.h */
+#ifndef XSRETURN_UV /* Er, almost. Fix 21086 from Sep 2003 */
+ #define XST_mUV(i,v) (ST(i) = sv_2mortal(newSVuv(v)) )
+ #define XSRETURN_UV(v) STMT_START { XST_mUV(0,v); XSRETURN(1); } STMT_END
+#endif
#include "ptypes.h"
#include "cache.h"
#include "sieve.h"
@@ -436,3 +440,26 @@ _XS_moebius(IN UV lo, IN UV hi = 0)
IV
_XS_mertens(IN UV n)
+
+UV
+_XS_exp_mangoldt(IN UV n)
+ CODE:
+ if (n <= 1) XSRETURN_UV(1);
+ if ((n & (n-1)) == 0) XSRETURN_UV(2); /* Power of 2 */
+ if ((n & 1) == 0) XSRETURN_UV(1); /* Even number */
+ /* if (_XS_is_prime(n)) XSRETURN_UV(n); */
+ {
+ UV factors[MPU_MAX_FACTORS+1];
+ UV nfactors, i;
+ /* We could try a partial factor, e.g. looking for two small factors */
+ /* We could also check powers of primes searching for n */
+ nfactors = factor(n, factors);
+ for (i = 1; i < nfactors; i++) {
+ if (factors[i] != factors[0])
+ XSRETURN_UV(1);
+ }
+ XSRETURN_UV(factors[0]);
+ }
+ RETVAL = 1;
+ OUTPUT:
+ RETVAL
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 1905f0c..f63495f 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -58,6 +58,7 @@ sub _import_nobigint {
undef *moebius; *moebius = \&_XS_moebius;
undef *mertens; *mertens = \&_XS_mertens;
undef *euler_phi; *euler_phi = \&_XS_totient;
+ undef *exp_mangoldt; *exp_mangoldt = \&_XS_exp_mangoldt;
}
BEGIN {
@@ -1172,7 +1173,7 @@ sub euler_phi {
next unless $totients[$i] == $i;
$totients[$i] = $i-1;
foreach my $j (2 .. int($nend / $i)) {
- $totients[$i*$j] = ($totients[$i*$j]*($i-1))/$i;
+ $totients[$i*$j] -= $totients[$i*$j]/$i;
}
}
splice(@totients, 0, $n) if $n > 0;
@@ -1180,28 +1181,30 @@ sub euler_phi {
}
return $n if $n <= 1;
- my %factor_mult;
- my @factors = grep { !$factor_mult{$_}++ } factor($n);
+ my @factors = factor($n);
if (ref($n) ne 'Math::BigInt') {
my $totient = 1;
- foreach my $factor (@factors) {
- $totient *= ($factor - 1);
- $totient *= $factor for (2 .. $factor_mult{$factor});
+ my $lastf = 0;
+ foreach my $f (@factors) {
+ if ($f == $lastf) { $totient *= $f; }
+ else { $totient *= $f-1; $lastf = $f; }
}
return $totient;
}
- # Some real wackiness to solve issues with Math::BigInt::GMP (not seen with
- # Pari or Calc). Results of the multiply will go negative if we don't do
- # this. To see if you hit the standalone bug:
- # perl -E 'my $a = 2931542417; use bigint lib=>'GMP'; my $n = 49754396241690624; my $x = $n*$a; say $x;'
- # This may be related to RT 71548 of Math::BigInt::GMP.
my $totient = $n->copy->bone;
+ my $lastf = 0;
foreach my $factor (@factors) {
+ # This screwball line is here to solve some issues with the GMP backend,
+ # which has a weird bug. Results of the multiply can turn negative (!)
+ # if we don't do this. Perhaps related to RT 71548?
+ # perl -le 'use Math::BigInt lib=>'GMP'; my $a = 2931542417; my $n = Math::BigInt->new("49754396241690624"); my $x = $n*$a; print $x;'
+ # perl -le 'use Math::BigInt lib=>'GMP'; my $a = Math::BigInt->bone; $a *= 2931542417; $a *= 49754396241690624; print $a;'
+ # TODO: more work reproducing this
my $f = $n->copy->bzero->badd("$factor");
- $totient->bmul($f->copy->bsub(1));
- $totient->bmul($f) for (2 .. $factor_mult{$factor});
+ if ($f == $lastf) { $totient->bmul($f); }
+ else { $totient->bmul($f->copy->bsub(1)); $lastf = $f; }
}
return $totient;
}
@@ -1279,16 +1282,17 @@ sub exp_mangoldt {
my($n) = @_;
return 1 if defined $n && $n <= 1;
_validate_positive_integer($n);
+ return _XS_exp_mangoldt($n) if $n <= $_XS_MAXVAL;
- #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);
+ # Power of 2
+ return 2 if ($n & ($n-1)) == 0;
+ # Even numbers can't be a power of an odd prime
+ return 1 unless $n & 1;
- return 1 unless scalar @factors == 1;
- return $factors[0];
+ my @factors = ($n <= $_XS_MAXVAL) ? _XS_factor($n) : factor($n);
+ my $first = shift @factors;
+ return 1 if scalar grep { $_ != $first } @factors;
+ return $first;
}
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 1542bc3..8081ee0 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -3,7 +3,8 @@ use strict;
use warnings;
use Test::More;
-use Math::Prime::Util qw/moebius mertens euler_phi jordan_totient divisor_sum/;
+use Math::Prime::Util
+ qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt/;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
@@ -109,6 +110,29 @@ my %sigmak = (
3 => [1, 9, 28, 73, 126, 252, 344, 585, 757, 1134, 1332, 2044, 2198, 3096, 3528, 4681, 4914, 6813, 6860, 9198, 9632, 11988, 12168, 16380, 15751, 19782, 20440, 25112, 24390, 31752, 29792, 37449, 37296, 44226, 43344, 55261, 50654, 61740, 61544],
);
+my %mangoldt = (
+ 0 => 1,
+ 1 => 1,
+ 2 => 2,
+ 3 => 3,
+ 4 => 2,
+ 5 => 5,
+ 6 => 1,
+ 7 => 7,
+ 8 => 2,
+ 9 => 3,
+ 10 => 1,
+ 11 => 11,
+ 25 => 5,
+ 27 => 3,
+ 399981 => 1,
+ 399982 => 1,
+ 399983 => 399983,
+ 823543 => 7,
+ 83521 => 17,
+ 130321 => 19,
+);
+
plan tests => 0 + 1
+ 1 # Small Moebius
@@ -120,7 +144,8 @@ plan tests => 0 + 1
+ 2 # Dedekind psi calculated two ways
+ 1 # Calculate J5 two different ways
+ 2 * $use64 # Jordan totient example
- + scalar(keys %sigmak);
+ + scalar(keys %sigmak)
+ + scalar(keys %mangoldt);
ok(!eval { moebius(0); }, "moebius(0)");
@@ -199,3 +224,8 @@ while (my($k, $sigmaref) = each (%sigmak)) {
}
is_deeply( \@slist, $sigmaref, "Sum of divisors to the ${k}th power: Sigma_$k" );
}
+
+###### Exponential of von Mangoldt
+while (my($n, $em) = each (%mangoldt)) {
+ is( exp_mangoldt($n), $em, "exp_mangoldt($n) == $em" );
+}
--
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