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