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