#define PERL_NO_GET_CONTEXT 1 /* Define at top for more efficiency. */
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "multicall.h" /* only works in 5.6 and newer */
/* We're not using anything for which we need ppport.h */
#ifndef XSRETURN_UV /* Er, almost. Fix 21086 from Sep 2003 */
#define XST_mUV(i,v) (ST(i) = sv_2mortal(newSVuv(v)) )
#define XSRETURN_UV(v) STMT_START { XST_mUV(0,v); XSRETURN(1); } STMT_END
#endif
#include "ptypes.h"
#include "cache.h"
#include "sieve.h"
#include "util.h"
#include "factor.h"
#include "lehmer.h"
#include "aks.h"
/* Workaround perl 5.6 UVs */
#if PERL_REVISION <= 5 && PERL_VERSION <= 6 && BITS_PER_WORD == 64
/* This could be blown up with a wacky string, but it's just for 5.6 */
#define set_val_from_sv(val, sv) \
{ char*ptr = SvPV_nolen(sv); val = Strtoul(ptr, NULL, 10); }
#else
#define set_val_from_sv(val, sv) \
val = SvUV(sv)
#endif
/* multicall compatibility stuff */
#if PERL_VERSION < 7
# define USE_MULTICALL 0 /* Too much trouble to work around it */
#else
# define USE_MULTICALL 1
#endif
#if PERL_VERSION < 13 || (PERL_VERSION == 13 && PERL_SUBVERSION < 9)
# define PERL_HAS_BAD_MULTICALL_REFCOUNT
#endif
#ifndef CvISXSUB
# define CvISXSUB(cv) CvXSUB(cv)
#endif
/* Not right, but close */
#if !defined cxinc && ( (PERL_VERSION == 8 && PERL_SUBVERSION >= 2) || (PERL_VERSION == 10 && PERL_SUBVERSION <= 1) )
# define cxinc() Perl_cxinc(aTHX)
#endif
static int pbrent_factor_a1(UV n, UV *factors, UV maxrounds) {
return pbrent_factor(n, factors, maxrounds, 1);
}
/* Is this a pedantically valid integer?
* Croaks if undefined or invalid.
* Returns 0 if it is an object or a string too large for a UV.
* Returns 1 if it is good to process by XS.
*/
static int _validate_int(SV* n, int negok)
{
dTHX;
char* ptr;
STRLEN i, len;
UV val;
int isneg = 0;
if (!SvOK(n)) croak("Parameter must be defined");
/* aside: to detect bigint: if ( SvROK(n) && sv_isa(n, "Math::BigInt") ) */
if (SvROK(n)) return 0;
/* Perhaps SvPVbyte, or other UTF8 stuff? */
ptr = SvPV(n, len);
if (len == 0) croak("Parameter '' must be a positive integer");
for (i = 0; i < len; i++) {
if (!isDIGIT(ptr[i])) {
if (i == 0 && ptr[i] == '-' && negok)
isneg = 1;
else if (i == 0 && ptr[i] == '+')
/* Allowed */ ;
else
croak("Parameter '%s' must be a positive integer", ptr); /* TODO NULL */
}
}
if (isneg) return -1; /* It's a valid negative number */
set_val_from_sv(val, n);
if (val == UV_MAX) { /* Could be bigger than UV_MAX. Need to find out. */
char vstr[40];
sprintf(vstr, "%"UVuf, val);
/* Skip possible leading zeros */
while (len > 0 && (*ptr == '0' || *ptr == '+'))
{ ptr++; len--; }
for (i = 0; i < len; i++)
if (vstr[i] != ptr[i])
return 0;
}
return 1;
}
/* Call a Perl sub with one SV* in and a SV* return value.
* We use this to foist off big inputs onto Perl.
*/
static SV* _callsub(SV* arg, const char* name)
{
dTHX;
dSP; /* Local copy of stack pointer */
int count;
SV* v;
ENTER; /* Start wrapper */
SAVETMPS; /* Start (2) */
PUSHMARK(SP); /* Start args: note our SP */
XPUSHs(arg);
PUTBACK; /* End args: set global SP to ours */
count = call_pv(name, G_SCALAR); /* Call the sub */
SPAGAIN; /* refresh local stack pointer */
if (count != 1)
croak("callback sub should return one value");
v = newSVsv(POPs); /* Get the returned value */
FREETMPS; /* End wrapper */
LEAVE; /* End (2) */
return v;
}
#if BITS_PER_WORD == 64
static const UV _max_prime = UVCONST(18446744073709551557);
#else
static const UV _max_prime = UVCONST(4294967291);
#endif
MODULE = Math::Prime::Util PACKAGE = Math::Prime::Util
PROTOTYPES: ENABLE
void
_XS_set_verbose(IN int verbose)
int
_XS_get_verbose()
void
_XS_set_callgmp(IN int call_gmp)
int
_XS_get_callgmp()
void
prime_precalc(IN UV n)
void
prime_memfree()
void
_prime_memfreeall()
UV
_XS_prime_count(IN UV low, IN UV high = 0)
CODE:
if (high == 0) { /* Without a Perl layer in front of this, we'll have */
high = low; /* the pathological case of a-0 turning into 0-a. */
low = 0;
}
if (GIMME_V == G_VOID) {
prime_precalc(high);
RETVAL = 0;
} else {
RETVAL = _XS_prime_count(low, high);
}
OUTPUT:
RETVAL
UV
_XS_legendre_pi(IN UV n)
UV
_XS_meissel_pi(IN UV n)
UV
_XS_lehmer_pi(IN UV n)
UV
_XS_LMO_pi(IN UV n)
UV
_XS_nth_prime(IN UV n)
int
_XS_is_aks_prime(IN UV n)
UV
_get_prime_cache_size()
CODE:
RETVAL = get_prime_cache(0, 0);
OUTPUT:
RETVAL
int
_XS_prime_maxbits()
CODE:
RETVAL = BITS_PER_WORD;
OUTPUT:
RETVAL
SV*
sieve_primes(IN UV low, IN UV high)
PREINIT:
AV* av = newAV();
CODE:
if (low <= high) {
START_DO_FOR_EACH_PRIME(low, high) {
av_push(av,newSVuv(p));
} END_DO_FOR_EACH_PRIME
}
RETVAL = newRV_noinc( (SV*) av );
OUTPUT:
RETVAL
SV*
trial_primes(IN UV low, IN UV high)
PREINIT:
UV curprime;
AV* av = newAV();
CODE:
if (low <= high) {
if (low >= 2) low--; /* Make sure low gets included */
curprime = _XS_next_prime(low);
while (curprime <= high) {
av_push(av,newSVuv(curprime));
curprime = _XS_next_prime(curprime);
}
}
RETVAL = newRV_noinc( (SV*) av );
OUTPUT:
RETVAL
SV*
segment_primes(IN UV low, IN UV high);
PREINIT:
AV* av = newAV();
CODE:
/* Could rewrite using {start/next/end}_segment_primes functions */
if ((low <= 2) && (high >= 2)) { av_push(av, newSVuv( 2 )); }
if ((low <= 3) && (high >= 3)) { av_push(av, newSVuv( 3 )); }
if ((low <= 5) && (high >= 5)) { av_push(av, newSVuv( 5 )); }
if (low < 7) low = 7;
if (low <= high) {
/* Call the segment siever one or more times */
UV low_d, high_d, segment_size;
unsigned char* sieve = get_prime_segment(&segment_size);
if (sieve == 0)
croak("Could not get segment cache");
/* To protect vs. overflow, work entirely with d. */
low_d = low / 30;
high_d = high / 30;
{ /* Avoid recalculations of this */
UV endp = (high_d >= (UV_MAX/30)) ? UV_MAX-2 : 30*high_d+29;
prime_precalc(isqrt(endp) + 1 );
}
while ( low_d <= high_d ) {
UV seghigh_d = ((high_d - low_d) < segment_size)
? high_d
: (low_d + segment_size-1);
UV range_d = seghigh_d - low_d + 1;
UV seghigh = (seghigh_d == high_d) ? high : (seghigh_d*30+29);
UV segbase = low_d * 30;
/* printf(" startd = %"UVuf" endd = %"UVuf"\n", startd, endd); */
MPUassert( seghigh_d >= low_d, "segment_primes highd < lowd");
MPUassert( range_d <= segment_size, "segment_primes range > segment size");
/* Sieve from startd*30+1 to endd*30+29. */
if (sieve_segment(sieve, low_d, seghigh_d) == 0) {
release_prime_segment(sieve);
croak("Could not segment sieve from %"UVuf" to %"UVuf, segbase+1, seghigh);
}
START_DO_FOR_EACH_SIEVE_PRIME( sieve, low - segbase, seghigh - segbase )
av_push(av,newSVuv( segbase + p ));
END_DO_FOR_EACH_SIEVE_PRIME
low_d += range_d;
low = seghigh+2;
}
release_prime_segment(sieve);
}
RETVAL = newRV_noinc( (SV*) av );
OUTPUT:
RETVAL
SV*
erat_primes(IN UV low, IN UV high)
PREINIT:
unsigned char* sieve;
AV* av = newAV();
CODE:
if (low <= high) {
sieve = sieve_erat30(high);
if (sieve == 0) {
croak("Could not generate sieve for %"UVuf, high);
} else {
if ((low <= 2) && (high >= 2)) { av_push(av, newSVuv( 2 )); }
if ((low <= 3) && (high >= 3)) { av_push(av, newSVuv( 3 )); }
if ((low <= 5) && (high >= 5)) { av_push(av, newSVuv( 5 )); }
if (low < 7) { low = 7; }
START_DO_FOR_EACH_SIEVE_PRIME( sieve, low, high ) {
av_push(av,newSVuv(p));
} END_DO_FOR_EACH_SIEVE_PRIME
Safefree(sieve);
}
}
RETVAL = newRV_noinc( (SV*) av );
OUTPUT:
RETVAL
#define SIMPLE_FACTOR(func, n, arg1) \
if (n <= 1) { \
XPUSHs(sv_2mortal(newSVuv( n ))); \
} else { \
while ( (n% 2) == 0 ) { n /= 2; XPUSHs(sv_2mortal(newSVuv( 2 ))); } \
while ( (n% 3) == 0 ) { n /= 3; XPUSHs(sv_2mortal(newSVuv( 3 ))); } \
while ( (n% 5) == 0 ) { n /= 5; XPUSHs(sv_2mortal(newSVuv( 5 ))); } \
if (n == 1) { /* done */ } \
else if (_XS_is_prime(n)) { XPUSHs(sv_2mortal(newSVuv( n ))); } \
else { \
UV factors[MPU_MAX_FACTORS+1]; \
int i, nfactors; \
nfactors = func(n, factors, arg1); \
for (i = 0; i < nfactors; i++) { \
XPUSHs(sv_2mortal(newSVuv( factors[i] ))); \
} \
} \
}
#define SIMPLE_FACTOR_2ARG(func, n, arg1, arg2) \
/* Stupid MSVC won't bring its C compiler out of the 1980s. */ \
if (n <= 1) { \
XPUSHs(sv_2mortal(newSVuv( n ))); \
} else { \
while ( (n% 2) == 0 ) { n /= 2; XPUSHs(sv_2mortal(newSVuv( 2 ))); } \
while ( (n% 3) == 0 ) { n /= 3; XPUSHs(sv_2mortal(newSVuv( 3 ))); } \
while ( (n% 5) == 0 ) { n /= 5; XPUSHs(sv_2mortal(newSVuv( 5 ))); } \
if (n == 1) { /* done */ } \
else if (_XS_is_prime(n)) { XPUSHs(sv_2mortal(newSVuv( n ))); } \
else { \
UV factors[MPU_MAX_FACTORS+1]; \
int i, nfactors; \
nfactors = func(n, factors, arg1, arg2); \
for (i = 0; i < nfactors; i++) { \
XPUSHs(sv_2mortal(newSVuv( factors[i] ))); \
} \
} \
}
void
_XS_factor(IN UV n)
PREINIT:
UV factors[MPU_MAX_FACTORS+1];
int i, nfactors;
PPCODE:
nfactors = factor(n, factors);
EXTEND(SP, nfactors);
for (i = 0; i < nfactors; i++) {
PUSHs(sv_2mortal(newSVuv( factors[i] )));
}
void
trial_factor(IN UV n, IN UV maxfactor = 0)
PPCODE:
SIMPLE_FACTOR(trial_factor, n, maxfactor);
void
fermat_factor(IN UV n, IN UV maxrounds = 64*1024*1024)
PPCODE:
SIMPLE_FACTOR(fermat_factor, n, maxrounds);
void
holf_factor(IN UV n, IN UV maxrounds = 8*1024*1024)
PPCODE:
SIMPLE_FACTOR(holf_factor, n, maxrounds);
void
squfof_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
PPCODE:
SIMPLE_FACTOR(squfof_factor, n, maxrounds);
void
rsqufof_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
PPCODE:
SIMPLE_FACTOR(racing_squfof_factor, n, maxrounds);
void
pbrent_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
PPCODE:
SIMPLE_FACTOR(pbrent_factor_a1, n, maxrounds);
void
prho_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
PPCODE:
SIMPLE_FACTOR(prho_factor, n, maxrounds);
void
pminus1_factor(IN UV n, IN UV B1 = 1*1024*1024, IN UV B2 = 0)
PPCODE:
if (B2 == 0)
B2 = 10*B1;
SIMPLE_FACTOR_2ARG(pminus1_factor, n, B1, B2);
int
_XS_is_pseudoprime(IN UV n, IN UV a)
int
_XS_miller_rabin(IN UV n, ...)
PREINIT:
UV bases[64];
int prob_prime = 1;
int c = 1;
CODE:
if (items < 2)
croak("No bases given to miller_rabin");
if ( (n == 0) || (n == 1) ) XSRETURN_IV(0); /* 0 and 1 are composite */
if ( (n == 2) || (n == 3) ) XSRETURN_IV(1); /* 2 and 3 are prime */
if (( n % 2 ) == 0) XSRETURN_IV(0); /* MR works with odd n */
while (c < items) {
int b = 0;
while (c < items) {
bases[b++] = SvUV(ST(c));
c++;
if (b == 64) break;
}
prob_prime = _XS_miller_rabin(n, bases, b);
if (prob_prime != 1)
break;
}
RETVAL = prob_prime;
OUTPUT:
RETVAL
int
_XS_is_extra_strong_lucas_pseudoprime(IN UV n)
int
_XS_is_frobenius_underwood_pseudoprime(IN UV n)
int
_XS_is_prob_prime(IN UV n)
int
_XS_is_prime(IN UV n)
int
is_prime(IN SV* n)
ALIAS:
is_prob_prime = 1
PREINIT:
int status;
CODE:
status = _validate_int(n, 1);
RETVAL = 0;
if (status == -1) {
/* return 0 */
} else if (status == 1) {
UV val;
set_val_from_sv(val, n);
RETVAL = _XS_is_prime(val);
} else {
SV* result;
const char* sub = 0;
if (_XS_get_callgmp())
sub = (ix == 0) ? "Math::Prime::Util::GMP::is_prime"
: "Math::Prime::Util::GMP::is_prob_prime";
else
sub = (ix == 0) ? "Math::Prime::Util::_generic_is_prime"
: "Math::Prime::Util::_generic_is_prob_prime";
result = _callsub(ST(0), sub);
RETVAL = SvIV(result);
}
OUTPUT:
RETVAL
UV
_XS_next_prime(IN UV n)
UV
_XS_prev_prime(IN UV n)
void
next_prime(IN SV* n)
ALIAS:
prev_prime = 1
PPCODE:
if (_validate_int(n, 0)) {
UV val;
set_val_from_sv(val, n);
if ( (ix && val < 3) || (!ix && val >= _max_prime) ) XSRETURN_UV(0);
if (ix) XSRETURN_UV(_XS_prev_prime(val));
else XSRETURN_UV(_XS_next_prime(val));
} else {
SV* result = _callsub(ST(0), (ix == 0) ?
"Math::Prime::Util::_generic_next_prime" :
"Math::Prime::Util::_generic_prev_prime" );
XPUSHs(sv_2mortal(result));
}
double
_XS_ExponentialIntegral(double x)
double
_XS_LogarithmicIntegral(double x)
double
_XS_RiemannZeta(double x)
CODE:
RETVAL = (double) ld_riemann_zeta(x);
OUTPUT:
RETVAL
double
_XS_RiemannR(double x)
void
_XS_totient(IN UV lo, IN UV hi = 0)
PREINIT:
UV i;
PPCODE:
if (hi != lo && hi != 0) { /* Totients in a range, returns array */
UV* totients;
if (hi < lo) XSRETURN_EMPTY;
if (lo < 2) {
if (lo <= 0 ) XPUSHs(sv_2mortal(newSVuv(0)));
if (lo <= 1 && hi >= 1) XPUSHs(sv_2mortal(newSVuv(1)));
lo = 2;
}
if (hi >= lo) {
totients = _totient_range(lo, hi);
/* Extend the stack to handle how many items we'll return */
EXTEND(SP, hi-lo+1);
for (i = lo; i <= hi; i++)
PUSHs(sv_2mortal(newSVuv(totients[i-lo])));
Safefree(totients);
}
} else {
UV facs[MPU_MAX_FACTORS+1]; /* maximum number of factors is log2n */
UV nfacs, totient, lastf;
UV n = lo;
if (n <= 1) XSRETURN_UV(n);
nfacs = factor(n, facs);
totient = 1;
lastf = 0;
for (i = 0; i < nfacs; i++) {
UV f = facs[i];
if (f == lastf) { totient *= f; }
else { totient *= f-1; lastf = f; }
}
PUSHs(sv_2mortal(newSVuv(totient)));
}
void
_XS_moebius(IN UV lo, IN UV hi = 0)
PREINIT:
UV i;
PPCODE:
if (hi != lo && hi != 0) { /* mobius in a range */
signed char* mu = _moebius_range(lo, hi);
MPUassert( mu != 0, "_moebius_range returned 0" );
EXTEND(SP, hi-lo+1);
for (i = lo; i <= hi; i++)
PUSHs(sv_2mortal(newSViv(mu[i-lo])));
Safefree(mu);
} else {
UV factors[MPU_MAX_FACTORS+1];
UV nfactors, lastf;
UV n = lo;
if (n <= 1)
XSRETURN_IV(n);
if ( (n >= 25) && ( !(n%4) || !(n%9) || !(n%25) ) )
XSRETURN_IV(0);
nfactors = factor(n, factors);
lastf = 0;
for (i = 0; i < nfactors; i++) {
if (factors[i] == lastf)
XSRETURN_IV(0);
lastf = factors[i];
}
XSRETURN_IV( (nfactors % 2) ? -1 : 1 );
}
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);
}
RETVAL = factors[0];
}
OUTPUT:
RETVAL
double
_XS_chebyshev_theta(IN UV n)
double
_XS_chebyshev_psi(IN UV n)
UV
_XS_divisor_sum(IN UV n)
int
_validate_num(SV* n, ...)
CODE:
RETVAL = 0;
if (_validate_int(n, 0)) {
if (items > 1 && SvOK(ST(1))) {
UV val, min, max;
set_val_from_sv(val, n);
set_val_from_sv(min, ST(1));
if (val < min)
croak("Parameter '%"UVuf"' must be >= %"UVuf, val, min);
if (items > 2 && SvOK(ST(2))) {
set_val_from_sv(max, ST(2));
if (val > max)
croak("Parameter '%"UVuf"' must be <= %"UVuf, val, max);
MPUassert( items <= 3, "_validate_num takes at most 3 parameters");
}
}
RETVAL = 1;
}
OUTPUT:
RETVAL
void
forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
PROTOTYPE: &$;$
CODE:
{
#if !USE_MULTICALL
dSP;
PUSHMARK(SP);
XPUSHs(block); XPUSHs(svbeg); XPUSHs(svend);
PUTBACK;
(void) call_pv("Math::Prime::Util::_generic_forprimes", G_VOID|G_DISCARD);
SPAGAIN;
#else
UV beg, end;
GV *gv;
HV *stash;
CV *cv;
SV* svarg;
void* ctx;
unsigned char* segment;
UV seg_base, seg_low, seg_high;
if (!_validate_int(svbeg, 0) || (items >= 3 && !_validate_int(svend,0))) {
dSP;
PUSHMARK(SP);
XPUSHs(block); XPUSHs(svbeg); XPUSHs(svend);
PUTBACK;
(void) call_pv("Math::Prime::Util::_generic_forprimes", G_VOID|G_DISCARD);
SPAGAIN;
XSRETURN_UNDEF;
}
if (items < 3) {
beg = 2;
set_val_from_sv(end, svbeg);
} else {
set_val_from_sv(beg, svbeg);
set_val_from_sv(end, svend);
}
if (beg > end)
XSRETURN_UNDEF;
cv = sv_2cv(block, &stash, &gv, 0);
if (cv == Nullcv)
croak("Not a subroutine reference");
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
/* Handle early part */
while (beg < 7) {
dSP;
beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
if (beg <= end) {
sv_setuv(svarg, beg);
GvSV(PL_defgv) = svarg;
PUSHMARK(SP);
call_sv((SV*)cv, G_VOID|G_DISCARD);
}
beg += 1 + (beg > 2);
}
if (beg <= end) {
ctx = start_segment_primes(beg, end, &segment);
while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
/* I'm getting a memory leak in the MULTICALL and I'm having no luck
* finding out why it is happening. Don't use this for now. */
if (0 && !CvISXSUB(cv)) {
dMULTICALL;
I32 gimme = G_VOID;
PUSH_MULTICALL(cv);
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) {
sv_setuv(svarg, seg_base + p);
GvSV(PL_defgv) = svarg;
MULTICALL;
} END_DO_FOR_EACH_SIEVE_PRIME
#ifdef PERL_HAS_BAD_MULTICALL_REFCOUNT
if (CvDEPTH(multicall_cv) > 1)
SvREFCNT_inc(multicall_cv);
#endif
POP_MULTICALL;
} else {
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) {
dSP;
sv_setuv(svarg, seg_base + p);
GvSV(PL_defgv) = svarg;
PUSHMARK(SP);
call_sv((SV*)cv, G_VOID|G_DISCARD);
} END_DO_FOR_EACH_SIEVE_PRIME
}
}
end_segment_primes(ctx);
}
SvREFCNT_dec(svarg);
#endif
XSRETURN_UNDEF;
}