[libmath-prime-util-perl] 18/35: exp_mangoldt in XS by default
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:50:03 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list