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 8aea368b9e763da2630463f2fa5812ea21485d76 Author: Dana Jacobsen <d...@acm.org> Date: Sat Mar 2 01:16:33 2013 -0800 10x speedup for divisor_sum --- Changes | 2 ++ XS.xs | 3 +++ factor.c | 24 ++++++++++++++++++++++++ factor.h | 2 ++ lib/Math/Prime/Util.pm | 37 +++++++++++++++++++++++++++++++------ 5 files changed, 62 insertions(+), 6 deletions(-) diff --git a/Changes b/Changes index e52cf01..4d3c052 100644 --- a/Changes +++ b/Changes @@ -20,6 +20,8 @@ Revision history for Perl extension Math::Prime::Util. - Add the first and second Chebyshev functions (theta and psi). + - Divisor sum with no sub is ~10x faster. + 0.22 26 February 2013 - Move main factor loop out of xs and into factor.c. diff --git a/XS.xs b/XS.xs index 83fc0ce..490b880 100644 --- a/XS.xs +++ b/XS.xs @@ -469,3 +469,6 @@ _XS_chebyshev_theta(IN UV n) double _XS_chebyshev_psi(IN UV n) + +UV +_XS_divisor_sum(IN UV n) diff --git a/factor.c b/factor.c index b86b02b..8bc1357 100644 --- a/factor.c +++ b/factor.c @@ -424,6 +424,30 @@ int _XS_is_prob_prime(UV n) } +UV _XS_divisor_sum(UV n) +{ + UV factors[MPU_MAX_FACTORS+1]; + int nfac, i; + UV product = 1; + + if (n <= 1) return n; + nfac = factor(n, factors); + for (i = 0; i < nfac; i++) { + if (i+1 < nfac && factors[i] == factors[i+1]) { + UV fmult = factors[i]*factors[i]; + do { + fmult *= factors[i++]; + } while (i+1 < nfac && factors[i] == factors[i+1]); + product *= (fmult-1) / (factors[i]-1); + } else { + product *= factors[i]+1; + } + } + return product; +} + + + /* Knuth volume 2, algorithm C. * Very fast for small numbers, grows rapidly. diff --git a/factor.h b/factor.h index b5cf735..a13b52e 100644 --- a/factor.h +++ b/factor.h @@ -25,6 +25,8 @@ extern int pminus1_factor(UV n, UV *factors, UV B1, UV B2); extern int _XS_miller_rabin(UV n, const UV *bases, int nbases); extern int _XS_is_prob_prime(UV n); +extern UV _XS_divisor_sum(UV n); + static UV gcd_ui(UV x, UV y) { UV t; if (y < x) { t = x; x = y; y = t; } diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index cc7129e..e491903 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -1289,16 +1289,26 @@ sub jordan_totient { # Mathematica and Pari both have functions like this. sub divisor_sum { my($n, $sub) = @_; - return 0 if defined $n && $n < 1; + # I really need to get cracking on an XS validator. + #return _XS_divisor_sum($n) if !defined $sub && defined $n && $n <= $_XS_MAXVAL && $_Config{'nobigint'}; + return (0,1)[$n] if defined $n && $n <= 1; _validate_positive_integer($n); if (!defined $sub) { - return $n if $n == 1; - my $sum = 1 + $n; - foreach my $f ( all_factors($n) ) { - $sum += $f; + return _XS_divisor_sum($n) if $n <= $_XS_MAXVAL; + my $product = 1; + my @factors = factor($n); + while (@factors) { + if (@factors > 1 && $factors[0] == $factors[1]) { + my $fmult = $factors[0] * $factors[0]; + $fmult *= shift @factors while @factors > 1 && $factors[0] == $factors[1]; + $product *= ($fmult -1) / ($factors[0] - 1); + } else { + $product *= $factors[0]+1; + } + shift @factors; } - return $sum; + return $product; } croak "Second argument must be a code ref" unless ref($sub) eq 'CODE'; @@ -3191,6 +3201,21 @@ Print some primes above 64-bit range: # Similar code using Pari: # perl -MMath::Pari=:int,PARI,nextprime -E 'my $start = PARI "100000000000000000000"; my $end = $start+1000; my $p=nextprime($start); while ($p <= $end) { say $p; $p = nextprime($p+1); }' +Examining the η3(x) function of Planat and Solé (2011): + + sub nu3 { + my $n = shift; + my $phix = chebyshev_psi($n); + my $nu3 = 0; + foreach my $nu (1..3) { + $nu3 += (moebius($nu)/$nu)*LogarithmicIntegral($phix**(1/$nu)); + } + return $nu3; + } + say prime_count(1000000); + say prime_count_approx(1000000); + say nu3(1000000); + Project Euler, problem 3 (Largest prime factor): use Math::Prime::Util qw/factor/; -- 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