[libmath-prime-util-perl] 75/181: More XS changes, plus small optimizations from bulk88
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:08 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit 517233e310bf4209ce5b25176941c0cd19be715f
Author: Dana Jacobsen <dana at acm.org>
Date: Tue Dec 31 01:31:24 2013 -0800
More XS changes, plus small optimizations from bulk88
---
XS.xs | 117 ++++++++++++++++++++-----------------------------
aks.c | 6 ++-
lib/Math/Prime/Util.pm | 8 ++--
util.c | 35 +++++++--------
util.h | 2 +-
5 files changed, 73 insertions(+), 95 deletions(-)
diff --git a/XS.xs b/XS.xs
index 851142f..2d8e507 100644
--- a/XS.xs
+++ b/XS.xs
@@ -129,9 +129,11 @@ static int _validate_int(pTHX_ SV* n, int negok)
static int _vcallsubn(pTHX_ I32 flags, const char* name, int nargs)
{
dSP;
+ char fullname[80] = "Math::Prime::Util::";
+ strncat(fullname, name, 60);
PUSHMARK(SP-nargs);
PUTBACK;
- return call_pv(name, flags);
+ return call_pv(fullname, flags);
}
#define _vcallsub(func) (void)_vcallsubn(aTHX_ G_SCALAR, func, 1)
@@ -163,7 +165,7 @@ prime_memfree()
case 2: ret = _XS_get_verbose(); break;
case 3: ret = _XS_get_callgmp(); break;
case 4: ret = get_prime_cache(0,0); break;
- case 5: ret = BITS_PER_WORD; break;
+ default: ret = BITS_PER_WORD; break;
}
if (ix > 1) XSRETURN_UV(ret);
@@ -223,7 +225,6 @@ _XS_nth_prime(IN UV n)
UV
_XS_legendre_phi(IN UV x, IN UV a)
-
SV*
sieve_primes(IN UV low, IN UV high)
ALIAS:
@@ -356,9 +357,6 @@ trial_factor(IN UV n, ...)
}
int
-_XS_is_pseudoprime(IN UV n, IN UV a)
-
-int
_XS_miller_rabin(IN UV n, ...)
PREINIT:
UV bases[64];
@@ -396,15 +394,17 @@ _XS_lucas_sequence(IN UV n, IN IV P, IN IV Q, IN UV k)
XPUSHs(sv_2mortal(newSVuv( Qk )));
int
-_XS_is_prime(IN UV n)
+_XS_is_prime(IN UV n, ...)
ALIAS:
_XS_is_prob_prime = 1
_XS_is_lucas_pseudoprime = 2
_XS_is_strong_lucas_pseudoprime = 3
_XS_is_extra_strong_lucas_pseudoprime = 4
_XS_is_frobenius_underwood_pseudoprime = 5
- _XS_is_aks_prime = 6
- _XS_BPSW = 7
+ _XS_is_almost_extra_strong_lucas_pseudoprime = 6
+ _XS_is_pseudoprime = 7
+ _XS_is_aks_prime = 8
+ _XS_BPSW = 9
PREINIT:
int ret;
CODE:
@@ -415,19 +415,16 @@ _XS_is_prime(IN UV n)
case 3: ret = _XS_is_lucas_pseudoprime(n, 1); break;
case 4: ret = _XS_is_lucas_pseudoprime(n, 2); break;
case 5: ret = _XS_is_frobenius_underwood_pseudoprime(n); break;
- case 6: ret = _XS_is_aks_prime(n); break;
+ case 6: ret = _XS_is_almost_extra_strong_lucas_pseudoprime
+ ( n, (items == 1) ? 1 : my_svuv(ST(1)) ); break;
+ case 7: ret = _XS_is_pseudoprime(n, my_svuv(ST(1))); break;
+ case 8: ret = _XS_is_aks_prime(n); break;
default:ret = _XS_BPSW(n); break;
}
RETVAL = ret;
OUTPUT:
RETVAL
-int
-_XS_is_almost_extra_strong_lucas_pseudoprime(IN UV n, IN UV increment = 1)
- CODE:
- RETVAL = _XS_is_almost_extra_strong_lucas_pseudoprime(n, increment);
- OUTPUT:
- RETVAL
void
is_prime(IN SV* svn)
@@ -445,11 +442,11 @@ is_prime(IN SV* svn)
} else {
const char* sub = 0;
if (_XS_get_callgmp())
- sub = (ix == 0) ? "Math::Prime::Util::GMP::is_prime"
- : "Math::Prime::Util::GMP::is_prob_prime";
+ sub = (ix == 0) ? "GMP::is_prime"
+ : "GMP::is_prob_prime";
else
- sub = (ix == 0) ? "Math::Prime::Util::_generic_is_prime"
- : "Math::Prime::Util::_generic_is_prob_prime";
+ sub = (ix == 0) ? "_generic_is_prime"
+ : "_generic_is_prob_prime";
_vcallsub(sub);
return; /* skip implicit PUTBACK */
}
@@ -468,8 +465,8 @@ next_prime(IN SV* svn)
XSRETURN_UV(_XS_next_prime(n));
}
}
- _vcallsub((ix == 0) ? "Math::Prime::Util::_generic_next_prime" :
- "Math::Prime::Util::_generic_prev_prime" );
+ _vcallsub((ix == 0) ? "_generic_next_prime" :
+ "_generic_prev_prime" );
return; /* skip implicit PUTBACK */
void
@@ -489,7 +486,7 @@ factor(IN SV* svn)
}
}
} else {
- _vcallsubn(aTHX_ GIMME_V, "Math::Prime::Util::_generic_factor", 1);
+ _vcallsubn(aTHX_ GIMME_V, "_generic_factor", 1);
return; /* skip implicit PUTBACK */
}
@@ -505,7 +502,7 @@ divisor_sum(IN SV* svn, ...)
UV sigma = divisor_sum(n, k);
if (sigma != 0) XSRETURN_UV(sigma); /* sigma 0 means overflow */
}
- _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_divisor_sum",items);
+ _vcallsubn(aTHX_ G_SCALAR, "_generic_divisor_sum",items);
return; /* skip implicit PUTBACK */
@@ -521,7 +518,7 @@ znorder(IN SV* sva, IN SV* svn)
if (order == 0) XSRETURN_UNDEF;
XSRETURN_UV(order);
}
- _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_znorder", 2);
+ _vcallsubn(aTHX_ G_SCALAR, "_generic_znorder", 2);
return; /* skip implicit PUTBACK */
void
@@ -536,7 +533,7 @@ znprimroot(IN SV* svn)
if (r == 0 && n != 1) XSRETURN_EMPTY;
XSRETURN_UV(r);
}
- _vcallsub("Math::Prime::Util::_generic_znprimroot");
+ _vcallsub("_generic_znprimroot");
return; /* skip implicit PUTBACK */
void
@@ -554,7 +551,7 @@ kronecker(IN SV* sva, IN SV* svb)
IV b = my_sviv(svb);
XSRETURN_IV( kronecker_ss(a, b) );
}
- _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_kronecker", 2);
+ _vcallsubn(aTHX_ G_SCALAR, "_generic_kronecker", 2);
return; /* skip implicit PUTBACK */
double
@@ -591,7 +588,7 @@ euler_phi(IN SV* svlo, ...)
UV lo = my_svuv(svlo);
XSRETURN_UV(totient(lo));
} else {
- _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_euler_phi", 1);
+ _vcallsubn(aTHX_ G_SCALAR, "_generic_euler_phi", 1);
return; /* skip implicit PUTBACK */
}
} else if (items == 2) {
@@ -617,7 +614,7 @@ euler_phi(IN SV* svlo, ...)
Safefree(totients);
}
} else {
- _vcallsubn(aTHX_ G_ARRAY,"Math::Prime::Util::_generic_euler_phi",items);
+ _vcallsubn(aTHX_ G_ARRAY,"_generic_euler_phi",items);
return; /* skip implicit PUTBACK */
}
} else {
@@ -633,7 +630,7 @@ moebius(IN SV* svlo, ...)
UV n = my_svuv(svlo);
XSRETURN_IV(moebius(n));
} else {
- _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_moebius",1);
+ _vcallsubn(aTHX_ G_SCALAR, "_generic_moebius",1);
return; /* skip implicit PUTBACK */
}
} else if (items == 2) {
@@ -655,34 +652,32 @@ moebius(IN SV* svlo, ...)
Safefree(mu);
}
} else {
- _vcallsubn(aTHX_ G_ARRAY, "Math::Prime::Util::_generic_moebius",items);
+ _vcallsubn(aTHX_ G_ARRAY, "_generic_moebius",items);
return; /* skip implicit PUTBACK */
}
} else {
croak("Usage: moebius(n) or moebius(1o,hi)");
}
-IV
-_XS_mertens(IN UV n)
-
-
void
carmichael_lambda(IN SV* svn)
ALIAS:
- exp_mangoldt = 1
- PREINIT:
- int status;
+ mertens = 1
+ exp_mangoldt = 2
PPCODE:
- status = _validate_int(aTHX_ svn, (ix == 0) ? 0 : 1);
- if (status == -1) {
- XSRETURN_UV(1);
- } else if (status == 1) {
- UV n = my_svuv(svn);
- if (ix == 0) XSRETURN_UV(carmichael_lambda(n));
- else XSRETURN_UV(exp_mangoldt(n));
+ int status = _validate_int(aTHX_ svn, (ix > 1) ? 1 : 0);
+ switch (ix) {
+ case 0: if (status == 1) XSRETURN_UV(carmichael_lambda(my_svuv(svn)));
+ _vcallsub("_generic_carmichael_lambda");
+ break;
+ case 1: if (status == 1) XSRETURN_IV(mertens(my_svuv(svn)));
+ _vcallsub("_generic_mertens");
+ break;
+ default:if (status ==-1) XSRETURN_UV(1);
+ if (status == 1) XSRETURN_UV(exp_mangoldt(my_svuv(svn)));
+ _vcallsub("_generic_exp_mangoldt");
+ break;
}
- _vcallsub( (ix == 0) ? "Math::Prime::Util::_generic_carmichael_lambda"
- : "Math::Prime::Util::_generic_exp_mangoldt" );
return; /* skip implicit PUTBACK */
int
@@ -717,11 +712,7 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
PPCODE:
{
#if !USE_MULTICALL
- PUSHMARK(SP);
- XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
- PUTBACK;
- (void) call_pv("Math::Prime::Util::_generic_forprimes", G_VOID|G_DISCARD);
- SPAGAIN;
+ _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_forprimes",items);
#else
UV beg, end;
GV *gv;
@@ -733,11 +724,7 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
UV seg_base, seg_low, seg_high;
if (!_validate_int(aTHX_ svbeg, 0) || (items >= 3 && !_validate_int(aTHX_ svend,0))) {
- PUSHMARK(SP);
- XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
- PUTBACK;
- (void) call_pv("Math::Prime::Util::_generic_forprimes", G_VOID|G_DISCARD);
- SPAGAIN;
+ _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_forprimes",items);
return;
}
if (items < 3) {
@@ -828,11 +815,7 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
if (items <= 1) return;
if (!_validate_int(aTHX_ svbeg, 0) || (items >= 3 && !_validate_int(aTHX_ svend,0))) {
- PUSHMARK(SP);
- XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
- PUTBACK;
- (void) call_pv("Math::Prime::Util::_generic_forcomposites", G_VOID|G_DISCARD);
- SPAGAIN;
+ _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_forcomposites",items);
return;
}
@@ -868,9 +851,9 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
ctx = start_segment_primes(beg, nextprime, &segment);
while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) {
- cbeg = (prevprime+1 < beg) ? beg : prevprime+1;
+ cbeg = prevprime+1; if (cbeg < beg) cbeg = beg;
prevprime = seg_base + p;
- cend = (prevprime-1 > end) ? end : prevprime-1;
+ cend = prevprime-1; if (cend > end) cbeg = end;
for (c = cbeg; c <= cend; c++) {
sv_setuv(svarg, c); MULTICALL;
}
@@ -913,11 +896,7 @@ fordivisors (SV* block, IN SV* svn)
if (items <= 1) return;
if (!_validate_int(aTHX_ svn, 0)) {
- PUSHMARK(SP);
- XPUSHs(block); XPUSHs(svn);
- PUTBACK;
- (void) call_pv("Math::Prime::Util::_generic_fordivisors", G_VOID|G_DISCARD);
- SPAGAIN;
+ _vcallsubn(aTHX_ G_VOID|G_DISCARD, "_generic_fordivisors",2);
return;
}
diff --git a/aks.c b/aks.c
index 6f4f9d8..733cb11 100644
--- a/aks.c
+++ b/aks.c
@@ -43,7 +43,6 @@ static int is_perfect_power(UV n) {
UV b, last;
if ((n <= 3) || (n == UV_MAX)) return 0;
if ((n & (n-1)) == 0) return 1; /* powers of 2 */
- last = log2floor(n-1) + 1;
#if (BITS_PER_WORD == 32) || (DBL_DIG > 19)
if (1) {
#elif DBL_DIG == 10
@@ -54,13 +53,16 @@ static int is_perfect_power(UV n) {
if ( n < (UV) pow(10, DBL_DIG) ) {
#endif
/* Simple floating point method. Fast, but need enough mantissa. */
- b = isqrt(n); if (b*b == n) return 1; /* perfect square */
+ b = isqrt(n);
+ if (b*b == n) return 1; /* perfect square */
+ last = log2floor(n-1) + 1;
for (b = 3; b < last; b = _XS_next_prime(b)) {
UV root = (UV) (pow(n, 1.0 / (double)b) + 0.5);
if ( ((UV)(pow(root, b)+0.5)) == n) return 1;
}
} else {
/* Dietzfelbinger, algorithm 2.3.5 (without optimized exponential) */
+ last = log2floor(n-1) + 1;
for (b = 2; b <= last; b++) {
UV a = 1;
UV c = n;
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 28f5a19..44c5b4d 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -66,8 +66,6 @@ sub _import_nobigint {
undef *nth_prime; *nth_prime = \&_XS_nth_prime;
undef *is_pseudoprime; *is_pseudoprime = \&_XS_is_pseudoprime;
undef *is_strong_pseudoprime; *is_strong_pseudoprime = \&_XS_miller_rabin;
- undef *moebius; *moebius = \&_XS_moebius;
- undef *mertens; *mertens = \&_XS_mertens;
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.
@@ -105,6 +103,7 @@ BEGIN {
*exp_mangoldt = \&Math::Prime::Util::_generic_exp_mangoldt;
*euler_phi = \&Math::Prime::Util::_generic_euler_phi;
*moebius = \&Math::Prime::Util::_generic_moebius;
+ *mertens = \&Math::Prime::Util::_generic_mertens;
*carmichael_lambda = \&Math::Prime::Util::_generic_carmichael_lambda;
*kronecker = \&Math::Prime::Util::_generic_kronecker;
*divisor_sum = \&Math::Prime::Util::_generic_divisor_sum;
@@ -1275,10 +1274,9 @@ sub _generic_moebius {
}
# A002321 Mertens' function. mertens(n) = sum(moebius(1,n))
-sub mertens {
+sub _generic_mertens {
my($n) = @_;
_validate_num($n) || _validate_positive_integer($n);
- return _XS_mertens($n) if $n <= $_XS_MAXVAL;
# This is the most basic Deléglise and Rivat algorithm. u = n^1/2
# and no segmenting is done. Their algorithm uses u = n^1/3, breaks
# the summation into two parts, and calculates those in segments. Their
@@ -4343,7 +4341,7 @@ Project Euler, problem 21 (Amicable numbers):
Project Euler, problem 41 (Pandigital prime), brute force command line:
- perl -MMath::Prime::Util=:all -E 'my @p = grep { /1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} @{primes(1000000,9999999)}; say $p[-1];'
+ perl -MMath::Prime::Util=primes -MList::Util=first -E 'say first { /1/&&/2/&&/3/&&/4/&&/5/&&/6/&&/7/} reverse @{primes(1000000,9999999)};'
Project Euler, problem 47 (Distinct primes factors):
diff --git a/util.c b/util.c
index 4f14cae..78d1714 100644
--- a/util.c
+++ b/util.c
@@ -197,17 +197,15 @@ static const unsigned char prime_sieve30[] =
/* Return of 2 if n is prime, 0 if not. Do it fast. */
int _XS_is_prime(UV n)
{
+ if (n <= 10)
+ return (n == 2 || n == 3 || n == 5 || n == 7) ? 2 : 0;
+
if (n < UVCONST(2000000000)) {
- UV d, m;
- unsigned char mtab;
+ UV d = n/30;
+ UV m = n - d*30;
+ unsigned char mtab = masktab30[m]; /* Bitmask in mod30 wheel */
const unsigned char* sieve;
- int isprime = -1;
-
- if (n <= 10) return (n == 2 || n == 3 || n == 5 || n == 7) ? 2 : 0;
-
- d = n/30;
- m = n - d*30;
- mtab = masktab30[m]; /* Bitmask in mod30 wheel */
+ int isprime;
/* Return 0 if a multiple of 2, 3, or 5 */
if (mtab == 0)
@@ -220,6 +218,7 @@ int _XS_is_prime(UV n)
if (!(n%7) || !(n%11) || !(n%13)) return 0;
/* Check primary cache */
+ isprime = -1;
if (n <= get_prime_cache(0, &sieve))
isprime = 2*((sieve[d] & mtab) == 0);
release_prime_cache(sieve);
@@ -635,12 +634,12 @@ static const unsigned short primes_small[] =
/* The nth prime will be less or equal to this number */
static UV _XS_nth_prime_upper(UV n)
{
- double fn = (double) n;
- double flogn, flog2n, upper;
+ double fn, flogn, flog2n, upper;
if (n < NPRIMES_SMALL)
return primes_small[n];
+ fn = (double) n;
flogn = log(n);
flog2n = log(flogn); /* Note distinction between log_2(n) and log^2(n) */
@@ -846,7 +845,7 @@ UV* _totient_range(UV lo, UV hi) {
return totients;
}
-IV _XS_mertens(UV n) {
+IV mertens(UV n) {
/* See Deléglise and Rivat (1996) for O(n^2/3 log(log(n))^1/3) algorithm.
* This implementation uses their lemma 2.1 directly, so is ~ O(n).
* In serial it is quite a bit faster than segmented summation of mu
@@ -948,11 +947,11 @@ int kronecker_ss(IV a, IV b) {
}
UV totient(UV n) {
- UV i, facs[MPU_MAX_FACTORS+1];
- UV nfacs = factor(n, facs);
- UV totient = 1;
- UV lastf = 0;
+ UV i, nfacs, totient, lastf, facs[MPU_MAX_FACTORS+1];
if (n <= 1) return n;
+ nfacs = factor(n, facs);
+ totient = 1;
+ lastf = 0;
for (i = 0; i < nfacs; i++) {
UV f = facs[i];
if (f == lastf) { totient *= f; }
@@ -1235,7 +1234,7 @@ double _XS_LogarithmicIntegral(double x) {
if (x == 0) return 0;
if (x == 1) return -INFINITY;
if (x == 2) return li2;
- if (x <= 0) croak("Invalid input to LogarithmicIntegral: x must be > 0");
+ if (x < 0) croak("Invalid input to LogarithmicIntegral: x must be >= 0");
return _XS_ExponentialIntegral(log(x));
}
@@ -1246,7 +1245,7 @@ UV _XS_Inverse_Li(UV x) {
UV lo = (UV) (n*logn);
UV hi = (UV) (n*logn * 2 + 2);
- if (x < 1) return 0;
+ if (x == 0) return 0;
if (hi <= lo) hi = UV_MAX;
while (lo < hi) {
UV mid = lo + (hi-lo)/2;
diff --git a/util.h b/util.h
index 6cb0cbb..0c40e59 100644
--- a/util.h
+++ b/util.h
@@ -17,7 +17,7 @@ extern UV _XS_nth_prime(UV x);
extern signed char* _moebius_range(UV low, UV high);
extern UV* _totient_range(UV low, UV high);
-extern IV _XS_mertens(UV n);
+extern IV mertens(UV n);
extern double _XS_chebyshev_theta(UV n);
extern double _XS_chebyshev_psi(UV n);
--
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