This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.33 in repository libmath-prime-util-perl.
commit 0ef064ed0fe54ea008aa4611e44025a6f4fe243b Author: Dana Jacobsen <d...@acm.org> Date: Tue Oct 29 20:51:39 2013 -0700 exp_mangoldt in XS by default --- Changes | 2 ++ XS.xs | 53 ++++++++++++++++++++++++------------------ lib/Math/Prime/Util.pm | 63 +++++++++++++++++++++++++------------------------- t/19-moebius.t | 1 + 4 files changed, 65 insertions(+), 54 deletions(-) diff --git a/Changes b/Changes index c923915..7e84105 100644 --- a/Changes +++ b/Changes @@ -19,6 +19,8 @@ Revision history for Perl module Math::Prime::Util - all_factors in scalar context returns sigma_0(n). + - exp_mangoldt defaults to XS for speed. + - Perl RiemannZeta changes: - Borwein Zeta calculations done in BigInt instead of BigFloat (speed). diff --git a/XS.xs b/XS.xs index c30c35b..f5b5f16 100644 --- a/XS.xs +++ b/XS.xs @@ -721,30 +721,39 @@ _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) - RETVAL = 1; - else if ((n & (n-1)) == 0) /* Power of 2 */ - RETVAL = 2; - else if ((n & 1) == 0) /* Even number */ - RETVAL = 1; - /* else if (_XS_is_prime(n)) RETVAL = n; */ - else { - 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); +void +exp_mangoldt(IN SV* svn) + PREINIT: + int status; + UV n; + PPCODE: + status = _validate_int(svn, 1); + if (status == -1) { + XSRETURN_UV(1); + } else if (status == 1) { + set_val_from_sv(n, svn); + if (n <= 1) + XSRETURN_UV(1); + else if ((n & (n-1)) == 0) /* Power of 2 */ + XSRETURN_UV(2); + else if ((n & 1) == 0) /* Even number */ + XSRETURN_UV(1); + else { + 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 = factors[0]; + } else { + _vcallsub("Math::Prime::Util::_generic_exp_mangoldt"); + XSRETURN(1); } - OUTPUT: - RETVAL int _validate_num(SV* n, ...) diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index ec194ea..3c8bca2 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -67,7 +67,6 @@ 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; undef *chebyshev_theta; *chebyshev_theta = \&_XS_chebyshev_theta; undef *chebyshev_psi; *chebyshev_psi = \&_XS_chebyshev_psi; # These should be fast anyway, but this skips validation. @@ -997,6 +996,7 @@ sub primes { # gets very slow as the bit size increases, but that is why we have the # method above for bigints. if (1) { + my $loop_limit = 2_000_000; if ($bits > $_Config{'maxbits'}) { my $p = Math::BigInt->bone->blsft($bits-1)->binc(); @@ -1012,29 +1012,32 @@ sub primes { } } croak "Random function broken?"; - } - # Send through the generic random_prime function. Decently fast, but - # quite a bit slower than the F&T A1 method above. - if (!defined $_random_nbit_ranges[$bits]) { - if ($bits > $_Config{'maxbits'}) { - my $low = Math::BigInt->new('2')->bpow($bits-1); - my $high = Math::BigInt->new('2')->bpow($bits); - # Don't pull the range in to primes, just odds - $_random_nbit_ranges[$bits] = [$low+1, $high-1]; - } else { - my $low = 1 << ($bits-1); - my $high = ($bits == $_Config{'maxbits'}) - ? ~0-1 - : ~0 >> ($_Config{'maxbits'} - $bits); - $_random_nbit_ranges[$bits] = [next_prime($low-1),prev_prime($high+1)]; - # Example: bits = 7. - # low = 1<<6 = 64. next_prime(64-1) = 67 - # high = ~0 >> (64-7) = 127. prev_prime(127+1) = 127 + } else { + + # Send through the generic random_prime function. Decently fast, but + # quite a bit slower than the F&T A1 method above. + if (!defined $_random_nbit_ranges[$bits]) { + if ($bits > $_Config{'maxbits'}) { + my $low = Math::BigInt->new('2')->bpow($bits-1); + my $high = Math::BigInt->new('2')->bpow($bits); + # Don't pull the range in to primes, just odds + $_random_nbit_ranges[$bits] = [$low+1, $high-1]; + } else { + my $low = 1 << ($bits-1); + my $high = ($bits == $_Config{'maxbits'}) + ? ~0-1 + : ~0 >> ($_Config{'maxbits'} - $bits); + $_random_nbit_ranges[$bits] = [next_prime($low-1),prev_prime($high+1)]; + # Example: bits = 7. + # low = 1<<6 = 64. next_prime(64-1) = 67 + # high = ~0 >> (64-7) = 127. prev_prime(127+1) = 127 + } } + my ($low, $high) = @{$_random_nbit_ranges[$bits]}; + return $_random_prime->($low, $high); + } - my ($low, $high) = @{$_random_nbit_ranges[$bits]}; - return $_random_prime->($low, $high); } sub random_maurer_prime { @@ -1637,16 +1640,13 @@ sub prime_iterator_object { # Exponential of Mangoldt function (A014963). # Return p if n = p^m [p prime, m >= 1], 1 otherwise. -sub exp_mangoldt { +sub _generic_exp_mangoldt { my($n) = @_; - return 1 if defined $n && $n <= 1; _validate_num($n) || _validate_positive_integer($n); - #return _XS_exp_mangoldt($n) if $n <= $_XS_MAXVAL; - # 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 if $n <= 1; # n <= 1 + return 2 if ($n & ($n-1)) == 0; # n power of 2 + return 1 unless $n & 1; # even n can't be p^m my @pe = ($n <= $_XS_MAXVAL) ? _XS_factor_exp($n) : factor_exp($n); return 1 if scalar @pe > 1; @@ -2823,7 +2823,6 @@ Some of the functions, including: moebius mertens euler_phi - exp_mangoldt chebyshev_theta chebyshev_psi is_prime @@ -2841,7 +2840,7 @@ will turn off bigint support for those functions. Those functions will then go directly to the XS versions, which will speed up very small inputs a B<lot>. This is useful if you're using the functions in a loop, but since the difference is less than a millisecond, it's really not important in general. The last -four functions have shortcuts by default so will only skip validation. +five functions have shortcuts by default so will only skip validation. If you are using bigints, here are some performance suggestions: @@ -3338,8 +3337,8 @@ than Pari 2.1.7's C<is_prime(n,1)> which is the default for L<Math::Pari>. =head2 prime_certificate - my @cert = prime_certificate($n); - say verify_prime(@cert) ? "proven prime" : "not prime"; + my $cert = prime_certificate($n); + say verify_prime($cert) ? "proven prime" : "not prime"; Given a positive integer C<n> as input, returns a primality certificate as a multi-line string. If we could not prove C<n> prime, an empty diff --git a/t/19-moebius.t b/t/19-moebius.t index 9db14d4..54a529e 100644 --- a/t/19-moebius.t +++ b/t/19-moebius.t @@ -116,6 +116,7 @@ my %sigmak = ( ); my %mangoldt = ( +-13 => 1, 0 => 1, 1 => 1, 2 => 2, -- 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