Math-Prime-Util
view release on metacpan or search on metacpan
return 1;
}
return 0;
}
/******************************************************************************/
/******************************************************************************/
MODULE = Math::Prime::Util PACKAGE = Math::Prime::Util
PROTOTYPES: ENABLE
BOOT:
{
int i;
HV * stash = gv_stashpv("Math::Prime::Util", TRUE);
newCONSTSUB(stash, "_ivsize", newSViv(IVSIZE));
newCONSTSUB(stash, "_uvsize", newSViv(UVSIZE));
newCONSTSUB(stash, "_uvbits", newSViv(UVSIZE * 8));
newCONSTSUB(stash, "_nvsize", newSViv(NVSIZE));
newCONSTSUB(stash, "_nvmantbits", newSViv(NVMANTBITS));
newCONSTSUB(stash, "_nvmantdigits", newSViv((IV)((NVMANTBITS+1) / 3.322)));
newCONSTSUB(stash, "_XS_prime_maxbits", newSViv(BITS_PER_WORD));
newCONSTSUB(stash, "_XS_has_uint64", newSViv(HAVE_UINT64));
newCONSTSUB(stash, "_XS_has_uint128", newSViv(HAVE_UINT128));
#if HAVE_FACTOR128
newCONSTSUB(stash, "_XS_factor_bits", newSViv(128));
#else
newCONSTSUB(stash, "_XS_factor_bits", newSViv(BITS_PER_WORD));
#endif
(void) integer_complexity(0); /* Initialize the shared cache lock. */
boot_register_custom_ops(aTHX);
{
MY_CXT_INIT;
MY_CXT.MPUroot = stash;
MY_CXT.MPUGMP = gv_stashpv("Math::Prime::Util::GMP", TRUE);
MY_CXT.MPUPP = gv_stashpv("Math::Prime::Util::PP", TRUE);
for (i = 0; i <= CINTS; i++) {
MY_CXT.const_int[i] = newSViv(i-1);
SvREADONLY_on(MY_CXT.const_int[i]);
}
New(0, MY_CXT.randcxt, csprng_context_size(), char);
csprng_init_seed(MY_CXT.randcxt);
MY_CXT.forcount = 0;
MY_CXT.forexit = 0;
MY_CXT.bigintname = NULL;
MY_CXT.bigintstash = NULL;
}
}
#if defined(USE_ITHREADS) && defined(MY_CXT_KEY)
void
CLONE(...)
PREINIT:
int i;
SV* bigintclass;
PPCODE:
{
MY_CXT_CLONE; /* possible declaration */
MY_CXT.MPUroot = gv_stashpv("Math::Prime::Util", TRUE);
MY_CXT.MPUGMP = gv_stashpv("Math::Prime::Util::GMP", TRUE);
MY_CXT.MPUPP = gv_stashpv("Math::Prime::Util::PP", TRUE);
/* These should be shared between threads, but that's dodgy. */
for (i = 0; i <= CINTS; i++) {
MY_CXT.const_int[i] = newSViv(i-1);
SvREADONLY_on(MY_CXT.const_int[i]);
}
/* Make a new CSPRNG context for this thread */
New(0, MY_CXT.randcxt, csprng_context_size(), char);
csprng_init_seed(MY_CXT.randcxt);
/* NOTE: There is no thread destroy, so these never get freed... */
MY_CXT.forcount = 0;
MY_CXT.forexit = 0;
bigintclass = get_sv("Math::Prime::Util::_BIGINT", 0);
xs_set_bigint_class(aTHX_ bigintclass);
}
return; /* skip implicit PUTBACK, returning @_ to caller, more efficient */
#endif
void
END(...)
PREINIT:
dMY_CXT;
int i;
PPCODE:
_prime_memfreeall();
MY_CXT.MPUroot = NULL;
MY_CXT.MPUGMP = NULL;
MY_CXT.MPUPP = NULL;
for (i = 0; i <= CINTS; i++) {
SV * const sv = MY_CXT.const_int[i];
MY_CXT.const_int[i] = NULL;
SvREFCNT_dec_NN(sv);
} /* stashes are owned by stash tree, no refcount on them in MY_CXT */
csprng_clear(MY_CXT.randcxt);
Safefree(MY_CXT.randcxt); MY_CXT.randcxt = 0;
MY_CXT.forcount = 0;
MY_CXT.forexit = 0;
MY_CXT.bigintname = NULL;
MY_CXT.bigintstash = NULL;
return; /* skip implicit PUTBACK, returning @_ to caller, more efficient */
void csrand(IN SV* seed = 0)
PREINIT:
dMY_CXT;
STRLEN size;
unsigned char* data;
unsigned char gmpseed[64];
volatile unsigned char* p;
uint32_t size32, n, i;
PPCODE:
if (items == 0 || !SvOK(seed)) {
csprng_init_seed(MY_CXT.randcxt);
if (_XS_get_callgmp() >= 42) {
n = (uint32_t) sizeof(gmpseed);
if (get_entropy_bytes(n, gmpseed) != n)
croak("Failed to get entropy bytes for GMP CSPRNG seed");
SEED_GMP_CSPRNG(n, gmpseed);
p = (volatile unsigned char*) gmpseed;
i = n;
while (i-- > 0) *p++ = 0;
}
} else if (_XS_get_secure()) {
croak("secure option set, manual seeding disabled");
} else {
data = (unsigned char*) SvPV(seed, size);
size32 = size > (STRLEN)UINT32_MAX ? UINT32_MAX : (uint32_t)size;
csprng_seed(MY_CXT.randcxt, size32, data);
SEED_GMP_CSPRNG(size32, data);
}
XSRETURN(0);
UV srand(IN UV seedval = 0)
PREINIT:
dMY_CXT;
unsigned char seed[8];
uint32_t seedbytes;
CODE:
if (_XS_get_secure())
croak("secure option set, manual seeding disabled");
if (items == 0) {
if (get_entropy_bytes(sizeof(UV),(unsigned char*)&seedval) != sizeof(UV))
croak("Failed to get entropy bytes for srand");
}
seedbytes = csprng_srand(MY_CXT.randcxt, seedval, seed);
SEED_GMP_CSPRNG(seedbytes, seed);
RETVAL = seedval;
OUTPUT:
RETVAL
void irand()
ALIAS:
irand32 = 1
irand64 = 2
PREINIT:
dMY_CXT;
PPCODE:
if (ix == 0 || ix == 1) {
XSRETURN_UV( irand32(MY_CXT.randcxt) );
} else {
#if BITS_PER_WORD == 64
XSRETURN_UV( irand64(MY_CXT.randcxt) );
#else
UV hi = irand32(MY_CXT.randcxt);
UV lo = irand32(MY_CXT.randcxt);
RETURN_UV_UV(hi, lo);
#endif
}
NV drand(NV m = 0.0)
ALIAS:
rand = 1
PREINIT:
dMY_CXT;
CODE:
PERL_UNUSED_VAR(ix);
RETVAL = drand64(MY_CXT.randcxt);
if (m != 0) RETVAL *= m;
OUTPUT:
RETVAL
void random_bytes(IN SV* svn)
ALIAS:
entropy_bytes = 1
PREINIT:
int nstatus;
UV n;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (nstatus != 1 || n > MAX_RANDOM_BYTES)
croak("%s: input must be an integer between 0 and %"UVuf, SUBNAME, MAX_RANDOM_BYTES);
if (n == 0)
RETURN_SV(sv_2mortal(newSVpvs("")));
{
dMY_CXT;
char* sptr;
SV* svret = sv_2mortal(newSV(n));
SvPOK_only(svret);
SvCUR_set(svret, n);
sptr = SvPVX(svret);
if (ix == 0) {
csprng_rand_bytes(MY_CXT.randcxt, n, (unsigned char*)sptr);
} else {
UV got = get_entropy_bytes(n, (unsigned char*)sptr);
if (got != n) croak("%s: requested %"UVuf" bytes, got %"UVuf, SUBNAME, n, got);
}
sptr[n] = '\0';
RETURN_SV(svret);
}
UV _is_csprng_well_seeded()
ALIAS:
_XS_get_verbose = 1
_XS_get_callgmp = 2
_XS_get_nobigint = 3
_XS_get_secure = 4
_XS_set_secure = 5
_get_forexit = 6
_start_for_loop = 7
_get_prime_cache_size = 8
CODE:
switch (ix) {
case 0: { dMY_CXT; RETVAL = is_csprng_well_seeded(MY_CXT.randcxt); } break;
case 1: RETVAL = _XS_get_verbose(); break;
case 2: RETVAL = _XS_get_callgmp(); break;
case 3: RETVAL = _XS_get_nobigint(); break;
case 4: RETVAL = _XS_get_secure(); break;
case 5: _XS_set_secure(); RETVAL = 1; break;
case 6: { dMY_CXT; RETVAL = MY_CXT.forexit; } break;
case 7: { dMY_CXT; MY_CXT.forcount++; RETVAL = MY_CXT.forexit; MY_CXT.forexit = 0; } break;
case 8:
default: RETVAL = get_prime_cache(0,0); break;
}
OUTPUT:
RETVAL
void
_XS_set_bigint_class(IN SV* sv)
CODE:
xs_set_bigint_class(aTHX_ sv);
bool _validate_integer(SV* svn)
ALIAS:
_validate_integer_nonneg = 1
_validate_integer_positive = 2
_validate_integer_abs = 3
PREINIT:
uint32_t mask;
CODE:
switch (ix) {
case 0: mask = IFLAG_ANY; break;
case 1: mask = IFLAG_NONNEG; break;
case 2: mask = IFLAG_POS; break;
case 3: mask = IFLAG_ABS; break;
default: croak("_validate_integer unknown flag value");
}
RETVAL = xs_validate_integer_inplace(aTHX_ svn, mask);
OUTPUT:
RETVAL
void _canonicalized_integer(SV* svn)
PPCODE:
RETURN_SV(xs_to_canonical(aTHX_ svn));
void _canonicalize_integers(SV* svr)
PREINIT:
SV *target, *out;
PPCODE:
if (!SvROK(svr))
croak("_canonicalize_integers: expected scalar or array reference");
target = SvRV(svr);
if (SvTYPE(target) == SVt_PVAV) {
xs_aref_to_canonical(aTHX_ svr, "_canonicalize_integers");
} else if (xs_is_sv_scalar_ref(svr)) {
out = xs_to_canonical(aTHX_ target);
if (out != target)
sv_setsv(target, out);
} else {
croak("_canonicalize_integers: expected scalar or array reference");
}
XSRETURN(0);
void prime_memfree()
PREINIT:
dMY_CXT;
PPCODE:
prime_memfree();
/* (void) _vcallgmpsubn(aTHX_ G_VOID|G_DISCARD, "_GMP_memfree", 0, 49); */
if (MY_CXT.MPUPP != NULL) DISPATCHPP_RETURN_VOID();
XSRETURN(0);
void prime_precalc(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) != 1)
croak("prime_precalc: n must fit in native unsigned integer");
prime_precalc(n);
XSRETURN(0);
void _XS_set_verbose(IN SV* svn)
ALIAS:
_XS_set_callgmp = 1
_XS_set_nobigint = 2
_end_for_loop = 3
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) != 1)
croak("%s: n must fit in native unsigned integer", SUBNAME);
switch (ix) {
case 0: _XS_set_verbose(n); break;
case 1: _XS_set_callgmp(n); break;
case 2: _XS_set_nobigint(n); break;
case 3:
default: { dMY_CXT; MY_CXT.forcount--; MY_CXT.forexit = n > 0; } break;
}
XSRETURN(0);
void prime_count(IN SV* svlo, IN SV* svhi = 0)
ALIAS:
semiprime_count = 1
twin_prime_count = 2
ramanujan_prime_count = 3
perfect_power_count = 4
prime_power_count = 5
lucky_count = 6
PREINIT:
UV lo = 0, hi, count = 0;
PPCODE:
if ((items == 1 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items == 2 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) && _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG))) {
if (lo <= hi) {
switch (ix) {
case 0: count = prime_count_range(lo, hi); break;
case 1: count = semiprime_count_range(lo, hi); break;
case 2: count = twin_prime_count_range(lo, hi); break;
case 3: count = ramanujan_prime_count_range(lo, hi); break;
case 4: count = perfect_power_count_range(lo, hi); break;
case 5: count = prime_power_count_range(lo, hi); break;
case 6: count = lucky_count_range(lo, hi); break;
}
}
XSRETURN_UV(count);
}
DISPATCHPP_RETURN();
void prime_count_upper(IN SV* svn)
ALIAS:
prime_count_lower = 1
prime_count_approx = 2
prime_power_count_upper = 3
prime_power_count_lower = 4
prime_power_count_approx = 5
perfect_power_count_upper = 6
perfect_power_count_lower = 7
perfect_power_count_approx = 8
ramanujan_prime_count_upper = 9
ramanujan_prime_count_lower = 10
ramanujan_prime_count_approx = 11
twin_prime_count_approx = 12
semiprime_count_approx = 13
lucky_count_upper = 14
lucky_count_lower = 15
lucky_count_approx = 16
PREINIT:
UV n, ret;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
switch (ix) {
case 0: ret = prime_count_upper(n); break;
case 1: ret = prime_count_lower(n); break;
case 2: ret = prime_count_approx(n); break;
case 3: ret = prime_power_count_upper(n); break;
case 4: ret = prime_power_count_lower(n); break;
case 5: ret = prime_power_count_approx(n); break;
case 6: ret = perfect_power_count_upper(n); break;
case 7: ret = perfect_power_count_lower(n); break;
case 8: ret = perfect_power_count_approx(n); break;
case 9: ret = ramanujan_prime_count_upper(n); break;
case 10: ret = ramanujan_prime_count_lower(n); break;
case 11: ret = ramanujan_prime_count_approx(n); break;
case 12: ret = twin_prime_count_approx(n); break;
case 13: ret = semiprime_count_approx(n); break;
case 14: ret = lucky_count_upper(n); break;
case 15: ret = lucky_count_lower(n); break;
case 16:
default: ret = lucky_count_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void sum_primes(IN SV* svlo, IN SV* svhi = 0)
PREINIT:
UV lo = 2, hi;
#if HAVE_FACTOR128 && HAVE_SUM_PRIMES128
uint64_t lo64 = 2, hi64;
#endif
PPCODE:
if ((items == 1 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items == 2 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) && _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG))) {
UV count = 0;
/* 32/64-bit, Legendre or table-accelerated sieving. */
if (sum_primes(lo, hi, &count))
XSRETURN_UV(count);
/* If that didn't work, try the 128-bit version if supported. */
#if HAVE_FACTOR128 && HAVE_SUM_PRIMES128
{
uint128_t sum128;
if (sum_primes128(lo, hi, &sum128))
RETURN_U128(sum128);
}
#endif
}
#if HAVE_FACTOR128 && HAVE_SUM_PRIMES128
else if ((items == 1 && xs_sv_to_uint64(aTHX_ &hi64, svlo)) ||
(items == 2 && xs_sv_to_uint64(aTHX_ &lo64, svlo) && xs_sv_to_uint64(aTHX_ &hi64, svhi))) {
uint128_t sum128;
if (sum_primes128(lo64, hi64, &sum128))
RETURN_U128(sum128);
}
#endif
DISPATCHPP_RETURN();
void random_prime(IN SV* svlo, IN SV* svhi = 0)
PREINIT:
UV lo = 2, hi, ret;
dMY_CXT;
PPCODE:
if ((items == 1 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items == 2 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) && _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG))) {
ret = random_prime(MY_CXT.randcxt,lo,hi);
if (ret) XSRETURN_UV(ret);
else XSRETURN_UNDEF;
}
DISPATCHPP_RETURN();
void print_primes(IN SV* svlo, IN SV* svhi = 0, IN SV* svfd = 0)
PREINIT:
UV lo, hi, fd;
int status = 1, ifd, tfd, close_status, print_status = 1, saved_errno = 0;
PPCODE:
status &= _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG);
if (items > 1) status &= _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG);
if (items > 2) status &= _validate_and_set(&fd, aTHX_ svfd, IFLAG_NONNEG);
if (items == 1)
{ hi = lo; lo = 2; }
if (status == 0)
DISPATCHPP_RETURN_VOID();
if (items == 3) {
if (fd > INT_MAX) croak("print_primes: fd out of range");
ifd = fd;
} else {
ifd = fileno(stdout);
}
tfd = PerlLIO_dup(ifd);
if (tfd < 0)
croak("print_primes: open fd %d failed: %s", ifd, Strerror(errno));
if (lo <= hi) {
print_status = print_primes(lo, hi, tfd);
if (!print_status) saved_errno = errno;
}
close_status = PerlLIO_close(tfd);
if (!print_status)
croak("print_primes write error: %s", Strerror(saved_errno));
if (close_status != 0)
croak("print_primes: close fd %d failed: %s", ifd, Strerror(errno));
XSRETURN(0);
UV
_LMO_pi(IN UV n)
ALIAS:
_legendre_pi = 1
_meissel_pi = 2
_lehmer_pi = 3
_LMOS_pi = 4
_segment_pi = 5
PREINIT:
UV ret;
CODE:
switch (ix) {
case 0: ret = LMO_prime_count(n); break;
case 1: ret = legendre_prime_count(n); break;
case 2: ret = meissel_prime_count(n); break;
case 3: ret = lehmer_prime_count(n); break;
case 4: ret = LMOS_prime_count(n); break;
default:ret = segment_prime_count(2,n); break;
}
RETVAL = ret;
OUTPUT:
RETVAL
void
sieve_primes(IN UV low, IN UV high)
ALIAS:
trial_primes = 1
erat_primes = 2
segment_primes = 3
PREINIT:
AV* av;
PPCODE:
CREATE_RETURN_AV(av);
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) {
if (ix == 0) { /* Sieve with primary cache */
START_DO_FOR_EACH_PRIME(low, high) {
av_push(av,newSVuv(p));
} END_DO_FOR_EACH_PRIME
} else if (ix == 1) { /* Trial */
for (low = next_prime(low-1);
low <= high && low != 0;
low = next_prime(low) ) {
av_push(av,newSVuv(low));
}
} else if (ix == 2) { /* Erat with private memory */
unsigned char* sieve = sieve_erat30(high);
START_DO_FOR_EACH_SIEVE_PRIME( sieve, 0, low, high ) {
av_push(av,newSVuv(p));
} END_DO_FOR_EACH_SIEVE_PRIME
Safefree(sieve);
} else if (ix == 3) { /* Segment */
unsigned char* segment;
UV seg_base, seg_low, seg_high;
void* ctx = start_segment_primes(low, high, &segment);
while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_base, seg_low, seg_high )
av_push(av,newSVuv( p ));
END_DO_FOR_EACH_SIEVE_PRIME
}
end_segment_primes(ctx);
}
}
XSRETURN(1);
void primes(IN SV* svlo, IN SV* svhi = 0)
PREINIT:
AV* av;
UV lo = 0, hi, i;
int lostatus = 1, histatus = 0;
bool native_ok = FALSE;
PPCODE:
if (items == 1) {
histatus = _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG);
native_ok = histatus != 0;
} else if (items == 2) {
lostatus = _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG);
histatus = _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG);
native_ok = lostatus != 0 && histatus != 0;
}
if (native_ok) {
CREATE_RETURN_AV(av);
if ((lo <= 2) && (hi >= 2)) av_push(av, newSVuv( 2 ));
if ((lo <= 3) && (hi >= 3)) av_push(av, newSVuv( 3 ));
if ((lo <= 5) && (hi >= 5)) av_push(av, newSVuv( 5 ));
if (lo < 7) lo = 7;
if (lo <= hi) {
if ( hi-lo <= 10
|| (hi > 100000000UL && hi-lo <= 330)
|| (hi > 4000000000UL && hi-lo <= 1500)
) {
for (i = !(lo&1); i <= hi-lo; i += 2)
if (is_prime(lo+i))
av_push(av,newSVuv(lo+i));
} else if (hi < (65536*30) || hi <= get_prime_cache(0,0)) {
START_DO_FOR_EACH_PRIME(lo, hi) {
av_push(av,newSVuv(p));
} END_DO_FOR_EACH_PRIME
} else {
unsigned char* segment;
UV seg_base, seg_low, seg_high;
void* ctx = start_segment_primes(lo, hi, &segment);
while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
START_DO_FOR_EACH_SIEVE_PRIME(segment, seg_base, seg_low, seg_high)
av_push(av,newSVuv( p ));
END_DO_FOR_EACH_SIEVE_PRIME
}
end_segment_primes(ctx);
}
}
XSRETURN(1);
}
DISPATCHPP_RETURN();
void almost_primes(IN UV k, IN SV* svlo, IN SV* svhi = 0)
ALIAS:
omega_primes = 1
PREINIT:
AV* av;
UV lo = 1, hi, i, n, *S;
PPCODE:
if ((items == 2 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items >= 3 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) && _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG))) {
CREATE_RETURN_AV(av);
S = 0;
if (ix == 0) n = generate_almost_primes(&S, k, lo, hi);
else n = range_omega_prime_sieve(&S, k, lo, hi);
for (i = 0; i < n; i++)
av_push(av, newSVuv(S[i]));
if (S != 0) Safefree(S);
XSRETURN(1);
}
DISPATCHPP_RETURN();
void prime_powers(IN SV* svlo, IN SV* svhi = 0)
ALIAS:
twin_primes = 1
semi_primes = 2
ramanujan_primes = 3
PREINIT:
AV* av;
UV lo = 0, hi, i, num, *L;
PPCODE:
if ((items == 1 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items == 2 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) && _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG))) {
CREATE_RETURN_AV(av);
if (ix == 0) { /* Prime power */
if ((lo <= 2) && (hi >= 2)) av_push(av, newSVuv( 2 ));
if ((lo <= 3) && (hi >= 3)) av_push(av, newSVuv( 3 ));
if ((lo <= 4) && (hi >= 4)) av_push(av, newSVuv( 4 ));
if ((lo <= 5) && (hi >= 5)) av_push(av, newSVuv( 5 ));
} else if (ix == 1) { /* Twin */
if ((lo <= 3) && (hi >= 3)) av_push(av, newSVuv( 3 ));
if ((lo <= 5) && (hi >= 5)) av_push(av, newSVuv( 5 ));
} else if (ix == 2) { /* Semi */
if ((lo <= 4) && (hi >= 4)) av_push(av, newSVuv( 4 ));
if ((lo <= 6) && (hi >= 6)) av_push(av, newSVuv( 6 ));
} else if (ix == 3) { /* Ramanujan */
if ((lo <= 2) && (hi >= 2)) av_push(av, newSVuv( 2 ));
}
if (lo < 7) lo = 7;
if (lo <= hi) {
switch (ix) {
case 0: num = prime_power_sieve(&L,lo,hi); break;
case 1: num = range_twin_prime_sieve(&L,lo,hi); break;
case 2: num = range_semiprime_sieve(&L,lo,hi); break;
case 3: num = range_ramanujan_prime_sieve(&L,lo,hi); break;
default: num = 0; L = 0; break;
}
for (i = 0; i < num; i++)
av_push(av,newSVuv(L[i]));
Safefree(L);
}
XSRETURN(1);
}
DISPATCHPP_RETURN();
void
lucky_numbers(IN SV* svlo, IN SV* svhi = 0)
PREINIT:
AV* av;
UV lo = 0, hi, i, nlucky = 0;
PPCODE:
if ((items == 1 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items == 2 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) && _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG))) {
CREATE_RETURN_AV(av);
if (lo == 0 && hi <= UVCONST(4000000000)) {
uint32_t* lucky = lucky_sieve32(&nlucky, hi);
for (i = 0; i < nlucky; i++)
av_push(av,newSVuv(lucky[i]));
Safefree(lucky);
} else {
UV* lucky = lucky_sieve_range(&nlucky, lo, hi);
for (i = 0; i < nlucky; i++)
av_push(av,newSVuv(lucky[i]));
Safefree(lucky);
}
XSRETURN(1);
}
DISPATCHPP_RETURN();
void minimal_goldbach_pair(IN SV* svn)
ALIAS:
goldbach_pair_count = 1
PREINIT:
UV n, res;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
if (ix == 0) {
res = minimal_goldbach_pair(n);
if (res == 0) XSRETURN_UNDEF;
} else {
res = goldbach_pair_count(n);
}
XSRETURN_UV(res);
}
DISPATCHPP_RETURN();
void goldbach_pairs(IN SV* svn)
PREINIT:
size_t npairs, i;
UV n, *L;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) == 1) {
if (GIMME_V != G_ARRAY)
XSRETURN_UV(goldbach_pair_count(n));
L = goldbach_pairs(&npairs, n);
if (L == 0) XSRETURN_EMPTY;
EXTEND(SP, (EXTEND_TYPE)npairs);
for (i = 0; i < npairs; i++)
PUSHs(sv_2mortal(newSVuv(L[i])));
Safefree(L);
XSRETURN(npairs);
}
DISPATCHPP_RETURN();
void powerful_numbers(IN SV* svlo, IN SV* svhi = 0, IN SV* svk = 0)
PREINIT:
int kstatus = 1;
AV* av;
UV lo = 1, hi, i, k = 2, npowerful, *powerful;
PPCODE:
if (items >= 3)
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus == 1 &&
((items == 1 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items >= 2 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) &&
_validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG)))) {
CREATE_RETURN_AV(av);
powerful = powerful_numbers_range(&npowerful, lo, hi, k);
for (i = 0; i < npowerful; i++)
av_push(av,newSVuv(powerful[i]));
Safefree(powerful);
XSRETURN(1);
}
DISPATCHPP_RETURN();
void sieve_range(IN SV* svn, IN SV* svwidth, IN SV* svdepth)
PREINIT:
int status;
UV n, width, depth, lo, hi, *P, np, i;
PPCODE:
/* Return index of every n unless it is a composite with factor > depth */
if (_validate_and_set(&width, aTHX_ svwidth, IFLAG_NONNEG) != 1)
croak("sieve_range: width must fit in native unsigned integer");
if (_validate_and_set(&depth, aTHX_ svdepth, IFLAG_NONNEG) != 1)
croak("sieve_range: depth must fit in native unsigned integer");
status = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (width == 0)
RETURN_NOTHING();
if (status != 1 || n > UV_MAX-(width-1))
DISPATCHPP_RETURN();
lo = (n<2) ? 2 : n;
hi = n + width - 1;
np = range_partial_sieve(&P, lo, hi, depth);
if (GIMME_V != G_ARRAY) {
Safefree(P);
XSRETURN_UV(np);
}
EXTEND(SP, (EXTEND_TYPE)np);
for (i = 0; i < np; i++)
PUSHs(sv_2mortal(newSVuv(P[i] - n)));
Safefree(P);
void
sieve_prime_cluster(IN SV* svlo, IN SV* svhi, ...)
PREINIT:
uint32_t nc, cl[100];
UV i, lo, hi, cval, nprimes, *list;
int done, leading_zero;
PPCODE:
nc = 1;
leading_zero = 0;
if (items > 100) croak("sieve_prime_cluster: too many entries");
cl[0] = 0;
for (i = 2; i < (UV)items; i++) {
if (!_validate_and_set(&cval, aTHX_ ST(i), IFLAG_NONNEG))
croak("sieve_prime_cluster: cluster values must be standard integers");
if (i == 2 && cval == 0) { leading_zero = 1; continue; }
if (cval & 1) croak("sieve_prime_cluster: values must be even");
if (cval > 2147483647UL) croak("sieve_prime_cluster: values must be 31-bit");
if (cval <= cl[nc-1]) croak("sieve_prime_cluster: values must be increasing");
cl[nc++] = cval;
}
done = 0;
if (_validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) &&
_validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG)) {
list = sieve_cluster(lo, hi, nc, cl, &nprimes);
if (list != 0) {
done = 1;
if (GIMME_V != G_ARRAY) {
Safefree(list);
XSRETURN_UV(nprimes);
}
EXTEND(SP, (EXTEND_TYPE)nprimes);
for (i = 0; i < nprimes; i++)
PUSHs(sv_2mortal(newSVuv( list[i] )));
Safefree(list);
}
}
if (!done) {
/* PP removes an explicit zero before calling older GMP backends. */
DISPATCHPP_RETURN_GMPIF(!leading_zero);
}
void is_pseudoprime(IN SV* svn, ...)
ALIAS:
is_euler_pseudoprime = 1
is_strong_pseudoprime = 2
PREINIT:
int i, status, ret = 0;
UV n, base;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == 1) {
if (n < 3) {
ret = (n >= 2);
} else if (ix >= 1 && !(n&1)) {
ret = 0;
} else if (items == 1) {
ret = (ix == 0) ? is_pseudoprime(n, 2) :
(ix == 1) ? is_euler_pseudoprime(n, 2) :
is_strong_pseudoprime(n, 2);
} else {
for (i = 1, ret = 1; i < items && ret == 1; i++) {
status = _validate_and_set(&base, aTHX_ ST(i), IFLAG_NONNEG);
if (status != 1) break;
if (base < 2) croak("%s: invalid base: %"UVuf, SUBNAME, base);
ret = (ix == 0) ? is_pseudoprime(n, base) :
(ix == 1) ? is_euler_pseudoprime(n, base) :
is_strong_pseudoprime(n, base);
}
}
}
if (status != 0) RETURN_NPARITY(ret);
DISPATCHPP_RETURN();
void is_prime(IN SV* svn)
ALIAS:
is_prob_prime = 1
is_provable_prime = 2
is_bpsw_prime = 3
is_aks_prime = 4
is_lucas_pseudoprime = 5
is_strong_lucas_pseudoprime = 6
is_extra_strong_lucas_pseudoprime = 7
is_frobenius_underwood_pseudoprime = 8
is_frobenius_khashin_pseudoprime = 9
is_catalan_pseudoprime = 10
is_euler_plumb_pseudoprime = 11
is_ramanujan_prime = 12
is_semiprime = 13
is_chen_prime = 14
is_safe_prime = 15
is_mersenne_prime = 16
PREINIT:
int status, ret;
UV n;
PPCODE:
ret = 0;
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == 1) {
switch (ix) {
case 0: ret = 2*is_prime(n); break;
case 1: ret = 2*is_prob_prime(n); break;
case 2: ret = 2*is_prime(n); break;
case 3: ret = 2*BPSW(n); break;
case 4: ret = is_aks_prime(n); break;
case 5: ret = is_lucas_pseudoprime(n, 0); break;
case 6: ret = is_lucas_pseudoprime(n, 1); break;
case 7: ret = is_lucas_pseudoprime(n, 3); break;
case 8: ret = is_frobenius_underwood_pseudoprime(n); break;
case 9: ret = is_frobenius_khashin_pseudoprime(n); break;
case 10: ret = is_catalan_pseudoprime(n); break;
case 11: ret = is_euler_plumb_pseudoprime(n); break;
case 12: ret = is_ramanujan_prime(n); break;
case 13: ret = is_semiprime(n); break;
case 14: ret = is_chen_prime(n); break;
case 15: ret = is_safe_prime(n); break;
case 16: ret = is_mersenne_prime(n); if (ret == -1) status = 0; break;
default: break;
}
}
#if HAVE_FACTOR128
else {
int callgmp = _XS_get_callgmp();
uint128_t n128;
if ((ix == 1 || /* is_prob_prime => is_prime128 */
ix == 3 || /* is_bpsw_prime => is_bpsw128 */
(ix == 0 && callgmp < 1) || /* GMP tries harder, use if avail. */
(ix == 13 && callgmp < 42) /* GMP is faster for now. */
) && xs_sv_to_uint128(aTHX_ &n128, svn)) {
if (ix == 3) ret = is_bpsw128(n128);
else if (ix == 13) ret = is_semiprime128(n128);
else ret = is_prime128(n128);
status = 1;
}
}
#endif
if (status != 0) RETURN_NPARITY(ret);
DISPATCHPP_RETURN();
void
is_perrin_pseudoprime(IN SV* svn, IN SV* svk = 0)
ALIAS:
is_almost_extra_strong_lucas_pseudoprime = 1
is_delicate_prime = 2
PREINIT:
int nstatus, kstatus, ret = -1;
UV n, k;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items == 1) {
kstatus = 1;
k = (ix == 0) ? 0 : (ix == 1) ? 1 : 10;
} else {
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
}
if (kstatus == 1 && ix == 0) {
if (k > 3) croak("%s: restriction must be between 0 and 3", SUBNAME);
if (nstatus == 1) ret = is_perrin_pseudoprime(n, k);
}
if (kstatus == 1 && ix == 1) {
if (k < 1 || k > 256) croak("%s: invalid increment: %"UVuf, SUBNAME, k);
if (nstatus == 1) ret = is_almost_extra_strong_lucas_pseudoprime(n, k);
}
if (kstatus == 1 && ix == 2 && k <= UINT32_MAX) {
if (k < 2) croak("%s: invalid base: %"UVuf, SUBNAME, k);
if (nstatus == 1) ret = is_delicate_prime(n, (uint32_t)k);
}
if (kstatus == 1 && nstatus == -1) ret = 0; /* Negative n => 0 return */
if (ret >= 0)
RETURN_NPARITY(ret);
DISPATCHPP_RETURN_GMPIF(kstatus == 1);
void
is_frobenius_pseudoprime(IN SV* svn, IN SV* svp = 0, IN SV* svq = 0)
PREINIT:
int nstatus, pstatus, qstatus;
UV n;
IV P, Q, maxparam;
PPCODE:
if (items == 1) {
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (nstatus == -1) RETURN_NPARITY(0);
if (nstatus == 1) RETURN_NPARITY(is_frobenius_pseudoprime(n));
} else if (items == 3) {
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
pstatus = _validate_and_set((UV*)&P, aTHX_ svp, IFLAG_IV);
qstatus = _validate_and_set((UV*)&Q, aTHX_ svq, IFLAG_IV);
if (nstatus == -1) RETURN_NPARITY(0);
/* If |P| and |Q| are less than this, then D=P*P-4*Q cannot overflow IV */
maxparam = (BITS_PER_WORD == 64) ? (IV)UVCONST(3037000497) : 46338L;
if (nstatus != 0 && pstatus != 0 && qstatus != 0 &&
P <= maxparam && P >= -maxparam && Q <= maxparam && Q >= -maxparam)
RETURN_NPARITY(is_frobenius_pseudoprime_pq(n, P, Q));
} else
croak("is_frobenius_pseudoprime: expected 1 or 3 arguments");
DISPATCHPP_RETURN();
void
miller_rabin_random(IN SV* svn, IN SV* svnbases = 0)
PREINIT:
int nstatus, bstatus;
UV n, nbases;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items < 2) { bstatus = 1; nbases = 1; }
else bstatus=_validate_and_set(&nbases,aTHX_ svnbases,IFLAG_POS);
if (nstatus != 0 && bstatus != 0) {
dMY_CXT;
RETURN_NPARITY(nstatus == -1 ? 0 : is_mr_random(MY_CXT.randcxt,n,nbases));
}
DISPATCHPP_RETURN();
void is_gaussian_prime(IN SV* sva, IN SV* svb)
PREINIT:
UV a, b;
PPCODE:
if (_validate_and_set(&a, aTHX_ sva, IFLAG_ABS) &&
_validate_and_set(&b, aTHX_ svb, IFLAG_ABS)) {
if (a == 0) RETURN_NPARITY( ((b % 4) == 3) ? 2*is_prime(b) : 0 );
if (b == 0) RETURN_NPARITY( ((a % 4) == 3) ? 2*is_prime(a) : 0 );
if (a < HALF_WORD && b < HALF_WORD) {
UV aa = a*a, bb = b*b;
if (UV_MAX-aa >= bb)
RETURN_NPARITY( 2*is_prime(aa+bb) );
}
}
DISPATCHPP_RETURN();
void
gcd(...)
PROTOTYPE: @
ALIAS:
lcm = 1
PREINIT:
int i, status = 1;
UV ret, nullv, n;
PPCODE:
/* For each arg, while valid input, validate+gcd/lcm. Shortcut stop. */
if (ix == 0) { ret = 0; nullv = 1; }
else { ret = 1; nullv = 0; }
for (i = 0; i < items && ret != nullv && status != 0; i++) {
status = _validate_and_set(&n, aTHX_ ST(i), IFLAG_ABS);
if (status == 0) break;
if (i == 0) {
ret = n;
} else {
UV gcd = gcd_ui(ret, n);
if (ix == 0) {
ret = gcd;
} else {
n /= gcd;
if (n <= (UV_MAX / ret) ) ret *= n;
else status = 0; /* Overflow */
}
}
}
if (status != 0)
XSRETURN_UV(ret);
DISPATCHPP_RETURN();
void
vecmin(...)
PROTOTYPE: @
ALIAS:
vecmax = 1
PREINIT:
int i, status;
UV ret, n, retindex;
PPCODE:
if (items == 0) XSRETURN_UNDEF;
if (items == 1) RETURN_SV_CANONICAL(ST(0));
retindex = 0;
if ((status = _validate_and_set(&ret, aTHX_ ST(0), IFLAG_ANY)) != 0) {
int sign = status, minmax = (ix == 0);
for (i = 1; i < items; i++) {
status = _validate_and_set(&n, aTHX_ ST(i), IFLAG_ANY);
if (status == 0) break;
if (( (sign == -1 && status == 1) ||
(n >= ret && sign == status)
) ? !minmax : minmax ) {
sign = status;
ret = n;
retindex = i;
}
}
}
if (status != 0) {
/* retindex is already set. */
} else if (1) { /* Use string compares to decide the min/max */
int minmax = (ix == 0);
STRLEN alen, blen;
char *aptr, *bptr;
retindex = 0;
aptr = SvPV(ST(0), alen);
(void) strint_minmax(minmax, 0, 0, aptr, alen);
for (i = 1; i < items; i++) {
bptr = SvPV(ST(i), blen);
if (strint_minmax(minmax, aptr, alen, bptr, blen)) {
aptr = bptr;
alen = blen;
retindex = i;
}
}
} else {
DISPATCHPP_RETURN();
}
RETURN_SV_CANONICAL(ST(retindex));
void
vecsum(...)
PROTOTYPE: @
ALIAS:
vecprod = 1
PREINIT:
int i, status;
UV ret, n;
PPCODE:
if (items == 0)
XSRETURN_UV(ix == 0 ? 0 : 1);
status = 1;
if (ix == 0) {
UV lo = 0;
IV hi = 0;
for (ret = 0, i = 0; i < items; i++) {
status = _validate_and_set(&n, aTHX_ ST(i), IFLAG_ANY);
if (status == 0) break;
if (status == 1) hi += (n > (UV_MAX - lo));
else hi -= ((UV_MAX-n) >= lo);
lo += n;
}
if (status != 0)
RETURN_IV_UV(hi, lo);
} else if (ix == 1) {
int sign = 1;
for (ret = 1, i = 0; i < items; i++) {
status = _validate_and_set(&n, aTHX_ ST(i), IFLAG_ANY);
if (status == 0) break;
if (ret > 0 && n > UV_MAX/ret) { status = 0; break; }
sign *= status;
ret *= n;
}
if (status != 0 && sign == 1)
XSRETURN_UV(ret);
if (status != 0 && sign == -1 && ret <= (UV)IV_MAX)
XSRETURN_IV(neg_iv(ret));
}
if (_XS_get_callgmp() < 26) {
/* If we don't have GMP vecsum/vecprod, do it here. */
const char **sptr;
STRLEN *slen, rlen;
char *resstr;
New(0, sptr, items, const char*);
New(0, slen, items, STRLEN);
for (i = 0; i < items; i++) {
(void)_validate_int(aTHX_ ST(i), 1);
sptr[i] = SvPV_nomg(ST(i), slen[i]);
}
if (ix == 0) resstr = strint_vecsum( sptr, slen, items, &rlen);
else resstr = strint_vecprod(sptr, slen, items, &rlen);
Safefree(sptr);
Safefree(slen);
RETURN_SIGN_STRINT_STR(1, resstr, rlen);
}
DISPATCHPP_RETURN();
void
vecprefixsum(...)
PROTOTYPE: @
PREINIT:
int type;
size_t i, len;
UV *L;
PPCODE:
if (items == 0)
RETURN_NOTHING();
if (SvROK(ST(0)) && SvTYPE(SvRV(ST(0))) == SVt_PVAV) {
if (items != 1)
croak("vecprefixsum: expected integer list or single array reference");
type = arrayref_to_int_array(aTHX_ &len, &L, 0, ST(0), "vecprefixsum");
} else {
type = array_to_int_array(aTHX_ &len, &L, 0, &ST(0), items);
}
if (type == IARR_TYPE_NEG) {
IV *SL = (IV*)L;
for (i = 1; i < len; i++) {
IV a = SL[i-1], b = SL[i];
if (a > 0 && b > 0 && a > IV_MAX - b) break; /* overflow */
if (a < 0 && b < 0 && a < IV_MIN - b) break; /* underflow */
SL[i] = a + b;
}
if (i >= len) RETURN_LIST_VALS(len, L, 0);
} else if (type != IARR_TYPE_BAD) {
for (i = 1; i < len; i++) {
if (L[i] > UV_MAX - L[i-1])
break;
L[i] += L[i-1];
}
if (i >= len) RETURN_LIST_VALS(len, L, 1);
}
Safefree(L);
DISPATCHPP_RETURN();
void
vecextract(IN SV* x, IN SV* svm)
PREINIT:
AV* av;
UV mask, i = 0;
PPCODE:
CHECK_ARRAYREF(x);
av = (AV*) SvRV(x);
if (SvROK(svm) && SvTYPE(SvRV(svm)) == SVt_PVAV) {
SSize_t j;
DECL_ARREF(mav);
USE_ARREF(mav, svm, SUBNAME, AR_READ);
for (j = 0; (Size_t)j < len_mav; j++) {
int status = _validate_and_set(&mask, aTHX_ FETCH_ARREF(mav,j), IFLAG_IV);
REFRESH_ARREF(mav);
if (status != 0) {
SV** VV = av_fetch(av, (SSize_t)mask, 0);
XPUSHs( VV ? *VV : &PL_sv_undef );
} else {
croak("vecextract invalid index");
}
}
} else if (_validate_and_set(&mask, aTHX_ svm, IFLAG_NONNEG)) {
while (mask) {
if (mask & 1) {
SV **VV = av_fetch(av, i, 0);
XPUSHs( VV ? *VV : &PL_sv_undef );
}
i++;
mask >>= 1;
}
} else {
DISPATCHPP_RETURN();
}
void
vecequal(IN SV* a, IN SV* b)
PREINIT:
int res;
PPCODE:
res = _compare_array_refs(aTHX_ a, b);
if (res == AREF_CMP_DISPATCH)
DISPATCHPP_RETURN();
if (res == AREF_CMP_INVALID)
croak("vecequal: expected scalar or array reference");
RETURN_NPARITY(res);
void
vecmex(...)
ALIAS:
vecpmex = 1
PROTOTYPE: @
PREINIT:
char *setv;
int i, status = 1;
UV min, n;
uint32_t mask;
PPCODE:
if (ix == 0) {
min = 0;
mask = IFLAG_NONNEG;
} else {
min = 1;
mask = IFLAG_POS;
}
if (items == 0)
XSRETURN_UV(min);
Newz(0, setv, items, char);
for (i = 0; i < items; i++) {
status = _validate_and_set(&n, aTHX_ ST(i), mask);
/* Ignore any bigint */
if (status == 1 && n-min < (UV)items)
setv[n-min] = 1;
}
for (i = 0; i < items; i++)
if (setv[i] == 0)
break;
Safefree(setv);
XSRETURN_UV(i+min);
void
frobenius_number(...)
PROTOTYPE: @
PREINIT:
int i, found1 = 0;
UV fn, n, *A;
PPCODE:
if (items == 0) XSRETURN_UNDEF;
Newz(0, A, items, UV);
for (i = 0; i < items; i++) {
if (!_validate_and_set(&n, aTHX_ ST(i), IFLAG_POS)) break;
if (n == 1) found1 = 1;
A[i] = n;
}
if (i == items && !found1)
fn = frobenius_number(A, i);
Safefree(A);
if (i == items) {
if (found1) XSRETURN_IV(-1);
if (fn == 0) XSRETURN_UNDEF;
if (fn != UV_MAX) XSRETURN_UV(fn);
}
DISPATCHPP_RETURN();
void
chinese(...)
ALIAS:
chinese2 = 1
PROTOTYPE: @
PREINIT:
int i, status, astatus, nstatus;
UV ret, lcm, *an;
SV **psva, **psvn;
PPCODE:
status = 1;
New(0, an, 2*items, UV);
ret = 0;
for (i = 0; i < items; i++) {
AV* av;
CHECK_ARRAYREF(ST(i));
av = (AV*) SvRV(ST(i));
if (av_count(av) != 2) croak("%s: expected 2-element array reference",SUBNAME);
psva = av_fetch(av, 0, 0);
psvn = av_fetch(av, 1, 0);
if (psva == 0 || psvn == 0) { status = 0; break; }
astatus = _validate_and_set(an+i, aTHX_ *psva, IFLAG_ANY);
nstatus = _validate_and_set(an+i+items, aTHX_ *psvn, IFLAG_ABS);
if (astatus == 0 || nstatus == 0) { status = 0; break; }
if (an[i+items] == 0) {
XPUSHs(&PL_sv_undef);
if (ix == 1) XPUSHs(&PL_sv_undef);
XSRETURN(1 + ix);
}
_mod_with(an+i, astatus, an[i+items]);
}
if (status)
status = chinese(&ret, &lcm, an, an+items, items);
Safefree(an);
if (status) {
if (ix == 0) {
if (status < 0) XSRETURN_UNDEF;
else XSRETURN_UV(ret);
} else {
if (status < 0) {
XPUSHs(&PL_sv_undef);
XPUSHs(&PL_sv_undef);
} else {
XPUSHs(sv_2mortal(newSVuv( ret )));
XPUSHs(sv_2mortal(newSVuv( lcm )));
}
XSRETURN(2);
}
}
DISPATCHPP_RETURN();
void cornacchia(IN SV* svd, IN SV* svn)
PREINIT:
UV d, n, x, y;
PPCODE:
if (_validate_and_set(&d, aTHX_ svd, IFLAG_NONNEG) &&
_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) ) {
if (!cornacchia(&x, &y, d, n)) XSRETURN_UNDEF;
PUSHs(sv_2mortal(newSVuv( x )));
PUSHs(sv_2mortal(newSVuv( y )));
XSRETURN(2);
}
DISPATCHPP_RETURN();
void lucas_sequence(...)
PREINIT:
UV U, V, Qk, n, P, Q, k;
int nstatus, pstatus, qstatus, kstatus;
PPCODE:
if (items != 4) croak("lucas_sequence: n, P, Q, k");
nstatus = _validate_and_set(&n, aTHX_ ST(0), IFLAG_POS);
pstatus = _validate_and_set(&P, aTHX_ ST(1), IFLAG_IV);
qstatus = _validate_and_set(&Q, aTHX_ ST(2), IFLAG_IV);
kstatus = _validate_and_set(&k, aTHX_ ST(3), IFLAG_NONNEG);
if (nstatus && pstatus && qstatus && kstatus) {
lucas_seq(&U, &V, &Qk, n, (IV)P, (IV)Q, k);
PUSHs(sv_2mortal(newSVuv( U ))); /* 4 args in, 3 out, no EXTEND needed */
PUSHs(sv_2mortal(newSVuv( V )));
PUSHs(sv_2mortal(newSVuv( Qk )));
XSRETURN(3);
}
DISPATCHPP_RETURN();
void lucasuvmod(IN SV* svp, IN SV* svq, IN SV* svk, IN SV* svn)
ALIAS:
lucasumod = 1
lucasvmod = 2
PREINIT:
int pstatus, qstatus;
UV P, Q, k, n, U, V;
PPCODE:
pstatus = _validate_and_set(&P, aTHX_ svp, IFLAG_ANY);
qstatus = _validate_and_set(&Q, aTHX_ svq, IFLAG_ANY);
if ((pstatus != 0) && (qstatus != 0) &&
_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG) &&
_validate_and_set(&n, aTHX_ svn, IFLAG_ABS)
) {
if (n == 0) XSRETURN_UNDEF;
P = (pstatus == 1) ? P % n : ivmod((IV)P,n);
Q = (qstatus == 1) ? Q % n : ivmod((IV)Q,n);
if (ix == 1) XSRETURN_UV(lucasumod(P, Q, k, n));
if (ix == 2) XSRETURN_UV(lucasvmod(P, Q, k, n));
lucasuvmod(&U, &V, P, Q, k, n);
PUSHs(sv_2mortal(newSVuv( U )));
PUSHs(sv_2mortal(newSVuv( V )));
XSRETURN(2);
}
DISPATCHPP_RETURN();
void lucasuv(IN SV* svp, IN SV* svq, IN SV* svk)
ALIAS:
lucasu = 1
lucasv = 2
PREINIT:
UV k;
IV P, Q, U, V;
PPCODE:
if (_validate_and_set((UV*)&P, aTHX_ svp, IFLAG_IV) &&
_validate_and_set((UV*)&Q, aTHX_ svq, IFLAG_IV) &&
_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG) &&
lucasuv(&U, &V, P, Q, k)) {
if (ix == 1) XSRETURN_IV(U); /* U = lucasu(P,Q,k) */
if (ix == 2) XSRETURN_IV(V); /* V = lucasv(P,Q,k) */
PUSHs(sv_2mortal(newSViv( U ))); /* (U,V) = lucasuv(P,Q,k) */
PUSHs(sv_2mortal(newSViv( V )));
XSRETURN(2);
}
DISPATCHPP_RETURN();
void fibonacci(IN SV* svk)
ALIAS:
lucas_number = 1
PREINIT:
UV k, N;
int kstatus;
PPCODE:
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_ANY);
if (kstatus != 0) {
if (kstatus == -1) k = neg_iv(k);
N = ix == 0 ? fibonacci_number(k) : lucas_number(k);
if (k == 0 || N > 0) {
/* fib(-n) = -fib(n) for even n, luc(-n) = -luc(n) for odd n */
if (kstatus == 1 || k % 2 != (UV)ix)
XSRETURN_UV(N);
else if (N <= IV_MAX)
XSRETURN_IV(-(IV)N);
}
}
DISPATCHPP_RETURN();
void is_sum_of_squares(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus, ret;
UV n, k;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (items < 2) { kstatus = 1; k = 2; }
else { kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG); }
if (nstatus != 0 && kstatus != 0) {
switch (k) {
case 0: ret = (n==0); break;
case 1: ret = is_power(n,2); break;
case 2: ret = is_sum_of_two_squares(n); break;
case 3: ret = is_sum_of_three_squares(n); break;
default: ret = 1; break;
}
RETURN_NPARITY(ret);
}
DISPATCHPP_RETURN();
void is_square(IN SV* svn)
ALIAS:
is_carmichael = 1
is_quasi_carmichael = 2
is_perfect_power = 3
is_fundamental = 4
is_lucky = 5
is_practical = 6
is_perfect_number = 7
is_cyclic = 8
is_totient = 9
PREINIT:
int status, ret;
UV n;
PPCODE:
ret = 0;
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == 1) {
switch (ix) {
case 0: ret = is_perfect_square(n); break;
case 1: ret = is_carmichael(n); break;
case 2: ret = is_quasi_carmichael(n); break;
case 3: ret = is_perfect_power(n); break;
case 4: ret = is_fundamental(n,0); break;
case 5: ret = is_lucky(n); break;
case 6: ret = is_practical(n); break;
case 7: ret = is_perfect_number(n); break;
case 8: ret = is_cyclic(n); break;
case 9:
default:ret = is_totient(n); break;
}
} else if (status == -1) {
switch (ix) {
case 3: ret = is_perfect_power_neg(neg_iv(n)); break;
case 4: ret = is_fundamental(neg_iv(n),1); break;
default:break;
}
}
if (ix == 0 && status == 0 &&
xs_sv_is_perfect_square(aTHX_ svn, &ret))
RETURN_NPARITY(ret);
if (status != 0) RETURN_NPARITY(ret);
DISPATCHPP_RETURN();
void squarefree_kernel(IN SV* svn)
PREINIT:
int status;
UV n;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == -1)
XSRETURN_IV( neg_iv(squarefree_kernel(neg_iv(n))) );
if (status == 1)
XSRETURN_UV( squarefree_kernel(n) );
DISPATCHPP_RETURN();
void is_powerfree(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus = 1;
UV n, k = 2;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (items >= 2) kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus != 1 || k > UINT32_MAX) croak("%s: k must be <= 4294967295", SUBNAME);
if (nstatus != 0)
RETURN_NPARITY( is_powerfree(n,k) );
DISPATCHPP_RETURN();
void powerfree_count(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus = 1;
UV n, k = 2;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items >= 2) kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus != 1 || k > UINT32_MAX) croak("%s: k must be <= 4294967295", SUBNAME);
if (nstatus == -1)
XSRETURN_UV(0);
if (nstatus == 1)
XSRETURN_UV( powerfree_count(n,k) );
DISPATCHPP_RETURN();
void powerfree_sum(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus = 1;
UV n, k = 2, res;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items >= 2) kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus != 1 || k > UINT32_MAX) croak("%s: k must be <= 4294967295", SUBNAME);
if (nstatus == -1)
XSRETURN_UV(0);
if (nstatus == 1) {
res = powerfree_sum(n,k);
if (res != 0 || n == 0)
XSRETURN_UV(res);
/* res is 0 and n > 0, so we overflowed. Fall through to PP. */
}
DISPATCHPP_RETURN();
void powerfree_part(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus = 1;
UV n, k = 2;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items >= 2) kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus != 1 || k > UINT32_MAX) croak("%s: k must be <= 4294967295", SUBNAME);
if (nstatus != 0) {
if (nstatus == -1)
XSRETURN_IV( neg_iv(powerfree_part(neg_iv(n),k)) );
XSRETURN_UV( powerfree_part(n,k) );
}
DISPATCHPP_RETURN();
void powerfree_part_sum(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus = 1;
UV n, k = 2, res;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items >= 2) kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus != 1 || k > UINT32_MAX) croak("%s: k must be <= 4294967295", SUBNAME);
if (nstatus == -1)
XSRETURN_UV(0);
if (nstatus == 1) {
res = powerfree_part_sum(n,k);
if (res != 0 || n == 0)
XSRETURN_UV(res);
/* res is 0 and n > 0, so we overflowed. Fall through to PP. */
}
DISPATCHPP_RETURN();
void nth_powerfree(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus = 1;
UV n, k = 2, res;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (items >= 2) kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus != 1 || k > UINT32_MAX) croak("%s: k must be <= 4294967295", SUBNAME);
if (nstatus == 1) {
if (n == 0 || k < 2)
XSRETURN_UNDEF;
res = nth_powerfree(n,k);
if (res != 0)
XSRETURN_UV(res);
/* if res = 0, overflow */
}
DISPATCHPP_RETURN();
void
is_power(IN SV* svn, IN SV* svk = 0, IN SV* svroot = 0)
PREINIT:
int nstatus, kstatus, ret;
UV n, k;
uint32_t root;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items == 2 && xs_is_sv_scalar_ref(svk)) {
svroot = svk;
svk = 0;
}
if (svk) SvGETMAGIC(svk);
if (items < 2 || svk == 0 || !SvOK(svk)) {
kstatus = 1; k = 0;
} else {
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
}
if (items >= 3 && !xs_is_sv_scalar_ref(svroot))
croak("is_power: third argument not a scalar reference");
if (nstatus != 0 && kstatus != 0) {
if (k != 0) {
if (nstatus == -1) {
if (k % 2 == 0) RETURN_NPARITY(0); /* negative n even k return 0 */
n = neg_iv(n);
}
ret = is_power_ret(n, k, &root);
} else { /* k = 0 */
if (nstatus == -1)
n = neg_iv(n);
/* Following Pari/GP: ispower(0) = ispower(1) = ispower(-1) = 0 */
ret = (n <= 1) ? 0 : powerof_ret(n, &root);
if (nstatus == -1 && ret > 0 && ret % 2 == 0) {
uint32_t v = valuation(ret,2);
ret >>= v;
if (ret == 1) ret = 0;
if (ret) root = ipow(root,1U << v);
}
}
if (ret && svroot != 0) {
SV *svr = SvRV(svroot);
if (nstatus==1) sv_setuv(svr, k == 1 ? n : root);
else sv_setiv(svr, k == 1 ? (IV)neg_iv(n) : -(IV)root);
}
RETURN_NPARITY(ret);
}
DISPATCHPP_RETURN_GMPIF(svroot == 0);
void
is_prime_power(IN SV* svn, IN SV* svroot = 0)
PREINIT:
int status, ret;
UV n, root;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (items >= 2 && !xs_is_sv_scalar_ref(svroot))
croak("is_prime_power: second argument not a scalar reference");
if (status != 0) {
ret = (status == 1) ? prime_power(n, &root) : 0;
if (ret && items >= 2)
sv_setuv(SvRV(svroot), root);
RETURN_NPARITY(ret);
}
DISPATCHPP_RETURN_GMPIF(svroot == 0);
void
is_polygonal(IN SV* svn, IN UV k, IN SV* svroot = 0)
PREINIT:
UV n;
int status;
PPCODE:
if (k < 3) croak("is_polygonal: k must be >= 3");
if (items >= 3 && !xs_is_sv_scalar_ref(svroot))
croak("is_polygonal: third argument not a scalar reference");
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == -1)
RETURN_NPARITY(0);
if (status == 1) {
bool overflow = 0;
UV root = polygonal_root(n, k, &overflow);
UV result = (n == 0) || root;
if (!overflow) {
if (result && svroot != 0)
sv_setuv(SvRV(svroot), root);
RETURN_NPARITY(result);
}
}
DISPATCHPP_RETURN_GMPIF(svroot == 0);
void inverse_li(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
#if BITS_PER_WORD == 64
/* Larger values need more precision to identify the exact integer. */
if (n <= UVCONST(1000000000000) && n < MPU_MAX_PRIME_IDX)
#else
if (n < MPU_MAX_PRIME_IDX)
#endif
XSRETURN_UV(inverse_li(n));
}
DISPATCHPP_RETURN();
void inverse_li_nv(IN SV* svx)
PREINIT:
const char *xstr;
STRLEN xlen;
int xtype;
NV x, ret;
PPCODE:
SvGETMAGIC(svx);
if (SvOK(svx)) {
xstr = SvPV(svx, xlen);
xtype = grok_number(xstr, xlen, 0);
if (xtype && !(xtype & (IS_NUMBER_INFINITY | IS_NUMBER_NAN))) {
x = SvROK(svx) ? STRTONV(xstr) : SvNV(svx);
if (x >= 0.0 && MPU_NV_ISFINITE(x)) {
ret = (NV)ld_inverse_li(x,0);
if (MPU_NV_ISFINITE(ret))
XSRETURN_NV(ret);
croak("inverse_li_nv: result is outside native floating-point range");
}
}
}
croak("inverse_li_nv: x must be a finite non-negative real number");
void nth_prime(IN SV* svn)
ALIAS:
nth_prime_upper = 1
nth_prime_lower = 2
nth_prime_approx = 3
PREINIT:
UV n, ret;
PPCODE:
if ( _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
n <= MPU_MAX_PRIME_IDX ) {
if (n == 0) XSRETURN_UNDEF;
switch (ix) {
case 0: ret = nth_prime(n); break;
case 1: ret = nth_prime_upper(n); break;
case 2: ret = nth_prime_lower(n); break;
case 3:
default: ret = nth_prime_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void nth_prime_power(IN SV* svn)
ALIAS:
nth_prime_power_upper = 1
nth_prime_power_lower = 2
nth_prime_power_approx = 3
PREINIT:
UV n, ret;
PPCODE:
if ( _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
n <= MPU_MAX_PRIME_IDX ) {
if (n == 0) XSRETURN_UNDEF;
switch (ix) {
case 0: ret = nth_prime_power(n); break;
case 1: ret = nth_prime_power_upper(n); break;
case 2: ret = nth_prime_power_lower(n); break;
case 3:
default: ret = nth_prime_power_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void nth_perfect_power(IN SV* svn)
ALIAS:
nth_perfect_power_upper = 1
nth_perfect_power_lower = 2
nth_perfect_power_approx = 3
PREINIT:
UV n, ret;
PPCODE:
if ( _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
n <= MPU_MAX_PERFECT_POW_IDX ) {
if (n == 0) XSRETURN_UNDEF;
switch (ix) {
case 0: ret = nth_perfect_power(n); break;
case 1: ret = nth_perfect_power_upper(n); break;
case 2: ret = nth_perfect_power_lower(n); break;
case 3:
default: ret = nth_perfect_power_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void nth_ramanujan_prime(IN SV* svn)
ALIAS:
nth_ramanujan_prime_upper = 1
nth_ramanujan_prime_lower = 2
nth_ramanujan_prime_approx = 3
PREINIT:
UV n, ret;
PPCODE:
if ( _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
n <= MPU_MAX_RMJN_PRIME_IDX ) {
if (n == 0) XSRETURN_UNDEF;
switch (ix) {
case 0: ret = nth_ramanujan_prime(n); break;
case 1: ret = nth_ramanujan_prime_upper(n); break;
case 2: ret = nth_ramanujan_prime_lower(n); break;
case 3:
default: ret = nth_ramanujan_prime_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void nth_twin_prime(IN SV* svn)
ALIAS:
nth_twin_prime_approx = 1
PREINIT:
UV n, ret;
PPCODE:
if ( _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
n <= MPU_MAX_TWIN_PRIME_IDX ) {
if (n == 0) XSRETURN_UNDEF;
switch (ix) {
case 0: ret = nth_twin_prime(n); break;
case 1:
default: ret = nth_twin_prime_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void nth_semiprime(IN SV* svn)
ALIAS:
nth_semiprime_approx = 1
PREINIT:
UV n, ret;
PPCODE:
if ( _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
n <= MPU_MAX_SEMI_PRIME_IDX ) {
if (n == 0) XSRETURN_UNDEF;
switch (ix) {
case 0: ret = nth_semiprime(n); break;
case 1:
default: ret = nth_semiprime_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void nth_lucky(IN SV* svn)
ALIAS:
nth_lucky_upper = 1
nth_lucky_lower = 2
nth_lucky_approx = 3
PREINIT:
UV n, ret;
PPCODE:
if ( _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
n <= MPU_MAX_LUCKY_IDX ) {
if (n == 0) XSRETURN_UNDEF;
switch (ix) {
case 0: ret = nth_lucky(n); break;
case 1: ret = nth_lucky_upper(n); break;
case 2: ret = nth_lucky_lower(n); break;
case 3:
default: ret = nth_lucky_approx(n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void next_prime(IN SV* svn)
ALIAS:
prev_prime = 1
PREINIT:
UV n, ret;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)
&& !(ix == 0 && n >= MPU_MAX_PRIME)) {
ret = 0;
switch (ix) {
case 0: ret = next_prime(n); break;
case 1: ret = prev_prime(n); break;
default: break;
}
if (ret == 0) XSRETURN_UNDEF;
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void next_prime_power(IN SV* svn)
ALIAS:
prev_prime_power = 1
PREINIT:
UV n, ret;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_ABS)
&& !(ix == 0 && n >= MPU_MAX_PRIME)) {
ret = 0;
switch (ix) {
case 0: ret = next_prime_power(n); break;
case 1: ret = prev_prime_power(n); break;
default: break;
}
if (ret == 0) XSRETURN_UNDEF;
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void next_perfect_power(IN SV* svn)
PREINIT:
UV n;
int status;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == 1) {
n = next_perfect_power(n);
if (n != 0) XSRETURN_UV(n);
} else if (status == -1) { /* next perfect power: negative n */
n = next_perfect_power_neg(neg_iv(n));
XSRETURN_IV(neg_iv(n));
}
DISPATCHPP_RETURN();
void prev_perfect_power(IN SV* svn)
PREINIT:
UV n;
int status;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == 1) {
if (n == 0) XSRETURN_IV(-1);
n = prev_perfect_power(n);
XSRETURN_UV(n);
} else if (status == -1) { /* prev perfect power: negative n */
n = prev_perfect_power_neg(neg_iv(n));
if (n > 0 && n <= (UV)IV_MAX)
XSRETURN_IV(neg_iv(n));
}
DISPATCHPP_RETURN();
void next_chen_prime(IN SV* svn)
PREINIT:
UV n, ret;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
ret = next_chen_prime(n);
if (ret != 0) XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void urandomb(IN SV* svbits)
PREINIT:
UV bits;
int bstatus;
PPCODE:
bstatus = _validate_and_set(&bits, aTHX_ svbits, IFLAG_NONNEG);
if (bstatus != 1 || bits > MAX_RANDOM_BITS)
croak("%s: bits must be between 0 and %"UVuf, SUBNAME, MAX_RANDOM_BITS);
if (bits <= BITS_PER_WORD) {
dMY_CXT;
XSRETURN_UV( urandomb(MY_CXT.randcxt, bits) );
}
DISPATCHPP_RETURN();
void random_ndigit_prime(IN SV* svdigits)
PREINIT:
int dstatus;
UV digits;
PPCODE:
dstatus = _validate_and_set(&digits, aTHX_ svdigits, IFLAG_POS);
if (dstatus != 1 || digits < 1 || digits > MAX_RANDOM_DIGITS)
croak("%s: digits must be between 1 and %"UVuf, SUBNAME, MAX_RANDOM_DIGITS);
if (digits <= uvmax_maxlen) {
dMY_CXT;
UV res = random_ndigit_prime(MY_CXT.randcxt, digits);
if (res != 0) XSRETURN_UV(res);
}
/* We could implement the nobigint completely here if we cared */
DISPATCHPP_RETURN_GMPIF(!_XS_get_nobigint());
void random_nbit_prime(IN SV* svbits)
ALIAS:
random_shawe_taylor_prime = 1
random_maurer_prime = 2
random_proven_prime = 3
random_strong_prime = 4
random_safe_prime = 5
random_semiprime = 6
random_unrestricted_semiprime = 7
PREINIT:
UV bits, res;
static const unsigned int minbits[] = { 2, 2, 2, 2, 128, 3, 4, 3 };
int bstatus;
PPCODE:
bstatus = _validate_and_set(&bits, aTHX_ svbits, IFLAG_POS);
if (bstatus != 1 || bits < minbits[ix] || bits > MAX_RANDOM_BITS)
croak("%s: bits must be between %u and %"UVuf, SUBNAME, minbits[ix], MAX_RANDOM_BITS);
if (bits <= BITS_PER_WORD) {
dMY_CXT;
void *cs = MY_CXT.randcxt;
switch (ix) {
case 5: res = random_safe_prime(cs,bits); break;
case 6: res = random_semiprime(cs,bits); break;
case 7: res = random_unrestricted_semiprime(cs,bits); break;
default: res = random_nbit_prime(cs,bits); break;
}
if (res) XSRETURN_UV(res);
}
DISPATCHPP_RETURN();
void urandomm(IN SV* svn)
PREINIT:
UV n, ret;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_POS)) {
dMY_CXT;
ret = urandomm64(MY_CXT.randcxt, n);
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void urandomr(IN SV* svlo, IN SV* svhi)
PREINIT:
UV lo, hi;
int lo_s, hi_s;
PPCODE:
lo_s = _validate_and_set(&lo, aTHX_ svlo, IFLAG_ANY);
hi_s = _validate_and_set(&hi, aTHX_ svhi, IFLAG_ANY);
if (lo_s == 1 && hi_s == 1) {
dMY_CXT;
if (lo > hi) XSRETURN_UNDEF;
if (lo == hi) XSRETURN_UV(lo);
if (lo == 0 && hi == UV_MAX)
XSRETURN_UV( BITS_PER_WORD == 64 ? irand64(MY_CXT.randcxt)
: irand32(MY_CXT.randcxt) );
XSRETURN_UV(lo + urandomm64(MY_CXT.randcxt, hi-lo+1));
}
DISPATCHPP_RETURN();
void toint(IN SV* svn)
PREINIT:
const char *s;
STRLEN len;
uint32_t stype;
SV *numsv;
PPCODE:
SvGETMAGIC(svn);
if (!SvOK(svn)) XSRETURN_UV(0); /* undef returns 0 without warning */
/* Fastest path: if it is already a native int then we're done. */
/* TODO: verify 5.6.2 */
if (SVNUMTEST(svn)) {
if (SvIsUV(svn)) XSRETURN_UV(SvUVX(svn));
XSRETURN_IV(SvIVX(svn));
}
/* Simple NV within native UV/IV range: truncate toward zero. */
if (SvNOK(svn) && !SvROK(svn)) {
NV x = SvNV(svn);
/* This is a NV, there is no way to recover possible lost precision */
if (x >= 0.0 && x < (NV)UV_MAX + 1.0) XSRETURN_UV((UV)x);
if (x < 0.0 && x > (NV)IV_MIN - 1.0) XSRETURN_IV((IV)x);
}
/* Get the string and determine what kind of number it is */
s = SvPV(svn, len);
if (s == 0 || len == 0) XSRETURN_UV(0); /* Empty string return 0 */
numsv = _normalize_toint_string(aTHX_ svn, &s, &len);
stype = _parse_strnum(s, len);
if (stype & SNUMFLAG_INVALID)
croak("%s: '%" SVf "' is not a valid number", SUBNAME, svn);
if (stype == SNUMFLAG_NATIVE)
XSRETURN_UV(PSTRTOULL(s, NULL, 10));
if (stype == (SNUMFLAG_NEG|SNUMFLAG_NATIVE))
XSRETURN_IV(PSTRTOLL(s,NULL,10));
if (stype & SNUMFLAG_RADIX) {
UV n, base = (stype & SNUMFLAG_HEXSTR) ? 16 :
(stype & SNUMFLAG_OCTSTR) ? 8 : 2;
char *rstr = 0;
STRLEN rlen = 0;
int sign = (stype & SNUMFLAG_NEG) ? -1 : 1;
int status;
if (*s == '-' || *s == '+') { s++; len--; }
s += 2; len -= 2;
status = strint_radix_to_int(&n, &rstr, &rlen, s, len, base);
if (status == 1) {
if (sign > 0)
XSRETURN_UV(n);
RETURN_SIGNMAG_UV_UV(-1, 0, n);
}
if (status == 2)
RETURN_SIGN_STRINT_STR(sign, rstr, rlen);
croak("%s: internal radix conversion error", SUBNAME);
}
if (stype & SNUMFLAG_FP) {
NV x = SvNV(numsv); /* This might lose user precision, so check */
if (x >= 0.0 && x < nvuvmaxval+1.0) XSRETURN_UV((UV)x);
if (x < 0.0 && x > nvivminval-1.0) XSRETURN_IV((IV)x);
/* Doesn't fit, so go through Math::BigFloat */
RETURN_SV( xs_call_root_1_sv(aTHX_ "_int_from_float", numsv) );
}
if (stype & SNUMFLAG_BIGINT)
RETURN_SV(xs_to_canonical_bigint(aTHX_ numsv, s, len));
if (stype == SNUMFLAG_UNKNOWN)
RETURN_SV( xs_call_root_1_sv(aTHX_ "_int_from_float", numsv) );
croak("%s: internal numeric conversion error", SUBNAME);
void random_factored_integer(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_POS)) {
dMY_CXT;
int f, nf;
UV r, F[MPU_MAX_FACTORS+1];
AV* av = newAV();
r = random_factored_integer(MY_CXT.randcxt, n, &nf, F);
sort_uv_array(F, nf);
for (f = 0; f < nf; f++)
av_push(av, newSVuv(F[f]));
XPUSHs(sv_2mortal(newSVuv( r )));
XPUSHs(sv_2mortal(newRV_noinc( (SV*) av )));
} else {
DISPATCHPP_RETURN();
}
void contfrac(IN SV* svnum, IN SV* svden)
PREINIT:
UV num, den, *cf, rem;
int nstatus, dstatus, i, steps;
PPCODE:
nstatus = _validate_and_set(&num, aTHX_ svnum, IFLAG_ANY);
dstatus = _validate_and_set(&den, aTHX_ svden, IFLAG_POS);
if (nstatus == 0 || dstatus == 0)
DISPATCHPP_RETURN();
if (nstatus == -1) num = neg_iv(num);
steps = contfrac(&cf, &rem, num, den);
if (GIMME_V != G_ARRAY) {
int count = steps;
if (nstatus == -1 && steps > 1)
count += (cf[1] == 1) ? -1 : 1;
Safefree(cf);
XSRETURN_UV(count);
}
if (nstatus == 1) {
EXTEND(SP, (EXTEND_TYPE)steps);
for (i = 0; i < steps; i++)
PUSHs(sv_2mortal(newSVuv( cf[i] )));
} else {
IV negc0 = neg_iv(cf[0]);
EXTEND(SP, (EXTEND_TYPE)(steps + 1));
if (steps == 1) {
PUSHs(sv_2mortal(newSViv(negc0)));
} else if (cf[1] == 1) {
PUSHs(sv_2mortal(newSViv(negc0 - 1)));
PUSHs(sv_2mortal(newSVuv(cf[2] + 1)));
for (i = 3; i < steps; i++)
PUSHs(sv_2mortal(newSVuv( cf[i] )));
} else {
PUSHs(sv_2mortal(newSViv(negc0 - 1)));
PUSHs(sv_2mortal(newSVuv(1)));
PUSHs(sv_2mortal(newSVuv(cf[1] - 1)));
for (i = 2; i < steps; i++)
PUSHs(sv_2mortal(newSVuv( cf[i] )));
}
}
Safefree(cf);
void from_contfrac(...)
PROTOTYPE: @
PREINIT:
size_t i;
UV n, cfA0, cfA1, cfB0, cfB1, cfAn, cfBn;
int nstatus, overflow;
PPCODE:
nstatus = 1;
overflow = 0;
cfA0 = 1; cfA1 = 0;
cfB0 = 0; cfB1 = 1;
if (items > 0) {
nstatus = _validate_and_set(&n, aTHX_ ST(0), IFLAG_ANY);
/* TODO: handle negative n */
cfA1 = n;
for (i = 1; nstatus == 1 && i < (size_t) items; i++) {
if (!_validate_and_set(&n, aTHX_ ST(i), IFLAG_POS))
break;
/* check each step for overflow */
overflow = (UV_MAX/n < cfA1) || (UV_MAX/n < cfB1);
if (overflow) break;
cfAn = n * cfA1;
cfBn = n * cfB1;
overflow = (UV_MAX-cfAn < cfA0) || (UV_MAX-cfBn < cfB0);
if (overflow) break;
cfAn = cfAn + cfA0;
cfBn = cfBn + cfB0;
cfA0 = cfA1; cfA1 = cfAn;
cfB0 = cfB1; cfB1 = cfBn;
}
if (i < (size_t) items) /* Covers overflow */
nstatus = 0;
}
if (nstatus == 1) {
XPUSHs(sv_2mortal(newSVuv( cfA1 )));
XPUSHs(sv_2mortal(newSVuv( cfB1 )));
XSRETURN(2);
}
DISPATCHPP_RETURN();
void convergents(...)
PROTOTYPE: @
PREINIT:
size_t len;
UV *L, *P, *Q;
int type;
PPCODE:
if (items == 0)
RETURN_NOTHING();
type = array_to_int_array(aTHX_ &len, &L, 0, &ST(0), items);
/* We punt negative cases to PP */
if (!(type == IARR_TYPE_BAD || type == IARR_TYPE_NEG)) {
if (convergents(&P, &Q, L, len)) {
size_t i;
if (GIMME_V != G_ARRAY) {
Safefree(P);
Safefree(Q);
Safefree(L);
XSRETURN_UV(len);
}
EXTEND(SP, (EXTEND_TYPE)len);
for (i = 0; i < len; i++)
PUSH_2ELEM_AREF(P[i], Q[i]);
Safefree(P);
Safefree(Q);
Safefree(L);
XSRETURN(len);
}
}
Safefree(L);
DISPATCHPP_RETURN();
void bestrational(IN SV* svx, IN SV* svdbound)
PREINIT:
UV dbound, P, Q;
STRLEN xlen;
const char *xstr;
int xtype;
PPCODE:
SvGETMAGIC(svx);
if (_validate_and_set(&dbound, aTHX_ svdbound, IFLAG_POS) &&
SvOK(svx) &&
!sv_isobject(svx)) {
xstr = SvPV_nomg(svx, xlen);
xtype = grok_number(xstr, xlen, 0);
if (xtype & (IS_NUMBER_INFINITY | IS_NUMBER_NAN))
croak("bestrational: first argument must be a finite real number");
if (xtype) {
NV x = SvNV(svx);
if (MPU_NV_ISFINITE(x) && bestrational(&P, &Q, x, dbound)) {
if (x >= 0.0) {
XPUSHs(sv_2mortal(newSVuv(P)));
XPUSHs(sv_2mortal(newSVuv(Q)));
XSRETURN(2);
} else if (P <= IV_MAX) {
XPUSHs(sv_2mortal(newSViv(-(IV)P)));
XPUSHs(sv_2mortal(newSVuv(Q)));
XSRETURN(2);
}
}
}
}
DISPATCHPP_RETURN();
void next_calkin_wilf(IN SV* svnum, IN SV* svden)
ALIAS:
next_stern_brocot = 1
PREINIT:
UV num, den;
int status;
PPCODE:
if (_validate_and_set(&num, aTHX_ svnum, IFLAG_POS) && _validate_and_set(&den, aTHX_ svden, IFLAG_POS)) {
if (gcd_ui(num, den) != 1)
croak("%s: rational must be reduced", SUBNAME);
switch (ix) {
case 0: status = next_calkin_wilf(&num, &den); break;
case 1: status = next_stern_brocot(&num, &den); break;
default: status = 0; break;
}
if (status) {
XPUSHs(sv_2mortal(newSVuv( num )));
XPUSHs(sv_2mortal(newSVuv( den )));
XSRETURN(2);
}
}
DISPATCHPP_RETURN();
void calkin_wilf_n(IN SV* svnum, IN SV* svden)
ALIAS:
stern_brocot_n = 1
PREINIT:
UV num, den, n;
PPCODE:
if (_validate_and_set(&num, aTHX_ svnum, IFLAG_POS) && _validate_and_set(&den, aTHX_ svden, IFLAG_POS)) {
switch (ix) {
case 0: n = calkin_wilf_n(num, den); break;
case 1: n = stern_brocot_n(num, den); break;
default: n = 0; break;
}
if (n) XSRETURN_UV(n);
}
DISPATCHPP_RETURN();
void nth_calkin_wilf(IN SV* svn)
ALIAS:
nth_stern_brocot = 1
PREINIT:
UV n, num, den;
int status;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_POS)) {
switch (ix) {
case 0: status = nth_calkin_wilf(&num, &den, n); break;
case 1: status = nth_stern_brocot(&num, &den, n); break;
default: status = 0; break;
}
if (status) {
XPUSHs(sv_2mortal(newSVuv( num )));
XPUSHs(sv_2mortal(newSVuv( den )));
XSRETURN(2);
}
}
DISPATCHPP_RETURN();
void nth_stern_diatomic(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG))
XSRETURN_UV(nth_stern_diatomic(n));
DISPATCHPP_RETURN();
void farey(IN SV* svn, IN SV* svk = 0)
PREINIT:
UV n, k;
int wantsingle, kresult;
PPCODE:
wantsingle = svk != 0;
if (wantsingle) {
if (!_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG))
k = UV_MAX;
}
if (_validate_and_set(&n, aTHX_ svn, IFLAG_POS) && n <= UINT32_MAX) {
uint32_t n32 = (uint32_t) n;
if (wantsingle) {
uint32_t p, q;
kresult = kth_farey(n32, k, &p, &q);
if (kresult == 0) XSRETURN_UNDEF;
if (kresult == 1) {
PUSH_2ELEM_AREF(p, q);
XSRETURN(1);
}
} else if (GIMME_V != G_ARRAY) {
UV len = farey_length(n32);
if (len > 0)
XSRETURN_UV(len);
} else {
uint32_t *num, *den;
UV i, len = farey_array(n32, &num, &den);
if (len > 0) {
EXTEND(SP, (EXTEND_TYPE)len);
for (i = 0; i < len; i++)
PUSH_2ELEM_AREF(num[i], den[i]);
Safefree(num);
Safefree(den);
XSRETURN(len);
}
}
}
DISPATCHPP_RETURN();
void next_farey(IN SV* svn, IN SV* svfrac)
ALIAS:
farey_rank = 1
PREINIT:
SV **psvp, **psvq;
AV* av;
UV n, p64, q64;
int status;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_POS) && n <= UINT32_MAX) {
uint32_t n32 = (uint32_t) n;
CHECK_ARRAYREF(svfrac);
av = (AV*) SvRV(svfrac);
if (av_count(av) != 2) croak("%s: expected 2-element array reference", SUBNAME);
psvp = av_fetch(av, 0, 0);
psvq = av_fetch(av, 1, 0);
status = 1;
if (psvp == 0 || psvq == 0)
status = 0;
if (status != 0)
status = _validate_and_set(&p64, aTHX_ *psvp, IFLAG_NONNEG);
if (status != 0)
status = _validate_and_set(&q64, aTHX_ *psvq, IFLAG_POS);
if (status != 0 && p64 >= q64 && ix == 0)
XSRETURN_UNDEF;
if (status != 0 && p64 >= q64) {
UV len = farey_length(n32);
if (len > 0)
XSRETURN_UV(len - (p64 == q64));
status = 0; /* Force dispatch */
}
if (status != 0) {
UV g = gcd_ui(p64,q64);
if (g > 1) { p64 /= g; q64 /= g; }
if (p64 <= UINT32_MAX && q64 <= UINT32_MAX) {
uint32_t p = (uint32_t) p64, q = (uint32_t) q64;
if (ix == 1) {
UV rank = farey_rank(n32, p, q);
if (rank != UV_MAX)
XSRETURN_UV(rank);
} else {
if (next_farey(n32, &p, &q)) {
PUSH_2ELEM_AREF(p, q);
XSRETURN(1);
}
}
}
}
}
DISPATCHPP_RETURN();
void Pi(IN SV* svdigits = 0)
PREINIT:
#ifdef USE_QUADMATH
const UV mantsize = FLT128_DIG;
const NV pival = 3.141592653589793238462643383279502884197169Q;
#elif defined(USE_LONG_DOUBLE) && defined(HAS_LONG_DOUBLE)
const UV mantsize = LDBL_DIG;
const NV pival = 3.141592653589793238462643383279502884197169L;
#else
const UV mantsize = DBL_DIG;
const NV pival = 3.141592653589793238462643383279502884197169;
#endif
UV digits;
int status;
PPCODE:
status = 1;
digits = 0;
if (items > 0)
status = _validate_and_set(&digits, aTHX_ svdigits, IFLAG_NONNEG);
if (status != 1 || digits > mantsize)
DISPATCHPP_RETURN();
if (digits == 0)
XSRETURN_NV( pival );
{
char* out = pidigits(digits);
NV pi = STRTONV(out);
Safefree(out);
XSRETURN_NV( pi );
}
void bernfrac(IN SV* svn)
ALIAS:
harmfrac = 1
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) != 0) {
if (ix == 0) {
IV num; UV den;
if (bernfrac(&num, &den, n)) {
XPUSHs(sv_2mortal(newSViv( num )));
XPUSHs(sv_2mortal(newSVuv( den )));
XSRETURN(2);
}
} else {
UV num, den;
if (harmfrac(&num, &den, n)) {
XPUSHs(sv_2mortal(newSVuv( num )));
XPUSHs(sv_2mortal(newSVuv( den )));
XSRETURN(2);
}
}
}
DISPATCHPP_RETURN();
void
_pidigits(IN SV* svdigits)
PREINIT:
UV digits;
char* out;
PPCODE:
if (!_validate_and_set(&digits, aTHX_ svdigits, IFLAG_NONNEG) ||
digits > UINT32_MAX)
croak("_pidigits: input digits exceeds 32-bit limit");
if (digits == 0)
XSRETURN_PV("");
out = pidigits(digits);
XPUSHs(sv_2mortal(newSVpvn(out, digits + (digits>1))));
Safefree(out);
void inverse_totient(IN SV* svn)
PREINIT:
U32 gimme_v;
int status;
UV i, n, ntotients;
PPCODE:
gimme_v = GIMME_V;
status = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (status == 1) {
if (gimme_v == G_SCALAR) {
XSRETURN_UV( inverse_totient_count(n) );
} else if (gimme_v == G_ARRAY) {
UV* tots = inverse_totient_list(&ntotients, n);
if (ntotients != UV_MAX) {
EXTEND(SP, (EXTEND_TYPE)ntotients);
for (i = 0; i < ntotients; i++)
PUSHs(sv_2mortal(newSVuv( tots[i] )));
Safefree(tots);
XSRETURN(ntotients);
}
}
}
DISPATCHPP_RETURN();
void inverse_sigma0(IN SV* svk, IN SV* svlo, IN SV* svhi = 0)
ALIAS:
inverse_sigma0_count = 1
PREINIT:
AV* av;
UV k, lo = 1, hi, i, count, *list;
int kstatus, lostatus, histatus;
PPCODE:
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (items == 2) {
histatus = _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG);
lostatus = 1;
} else {
lostatus = _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG);
histatus = _validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG);
}
if (kstatus && lostatus && histatus) {
if (ix == 1)
XSRETURN_UV( inverse_sigma0_count(lo, hi, k) );
CREATE_RETURN_AV(av);
list = inverse_sigma0_list(&count, lo, hi, k);
for (i = 0; i < count; i++)
av_push(av, newSVuv(list[i]));
if (list != 0) Safefree(list);
XSRETURN(1);
}
DISPATCHPP_RETURN();
void
factor(IN SV* svn)
ALIAS:
factor_exp = 1
PREINIT:
UV n;
uint32_t i;
U32 gimme_v;
int status;
PPCODE:
gimme_v = GIMME_V;
status = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (status == 1) {
if (ix == 0) {
UV factors[MPU_MAX_FACTORS];
uint32_t nfactors = factor(n, factors);
if (gimme_v == G_SCALAR)
XSRETURN_UV(nfactors);
EXTEND(SP, (EXTEND_TYPE)nfactors);
for (i = 0; i < nfactors; i++)
PUSHs(sv_2mortal(newSVuv( factors[i] )));
} else {
factored_t nf = factorint(n);
if (gimme_v == G_SCALAR)
XSRETURN_UV(nf.nfactors);
EXTEND(SP, (EXTEND_TYPE)nf.nfactors);
for (i = 0; i < nf.nfactors; i++)
PUSH_2ELEM_AREF( nf.f[i], nf.e[i] );
}
} else {
#if HAVE_FACTOR128
if (_XS_get_callgmp() < 49) { /* Skip this if GMP backend will factor */
factored128_t nf;
if (xs_factorintp128_sv(aTHX_ &nf, svn)) {
uint32_t total;
uint16_t fi, ei;
if (ix == 0) {
/* flat list */
total = factored128p_total_factors(&nf);
if (gimme_v == G_SCALAR) XSRETURN_UV(total);
EXTEND(SP, (EXTEND_TYPE)total);
for (fi = 0; fi < nf.nfactors; fi++)
for (ei = 0; ei < nf.e[fi]; ei++)
PUSH_U128((uint128_t)nf.f[fi]);
if (nf.flarge)
PUSH_U128(nf.flarge);
} else {
/* [p, e] pairs */
total = factored128p_distinct_factors(&nf);
if (gimme_v == G_SCALAR) XSRETURN_UV(total);
EXTEND(SP, (EXTEND_TYPE)total);
for (fi = 0; fi < nf.nfactors; fi++) {
AV* av_ = newAV();
SV* sv_;
SV_FROM_U128(sv_, (uint128_t) nf.f[fi]);
av_push(av_, SvREFCNT_inc(sv_));
av_push(av_, newSVuv(nf.e[fi]));
PUSHs(sv_2mortal(newRV_noinc((SV*) av_)));
}
if (nf.flarge) {
AV* av_ = newAV();
SV* sv_;
SV_FROM_U128(sv_, nf.flarge);
av_push(av_, SvREFCNT_inc(sv_));
av_push(av_, newSVuv(1));
PUSHs(sv_2mortal(newRV_noinc((SV*) av_)));
}
}
XSRETURN(total);
}
}
#endif
DISPATCHPP_RETURN();
}
void divisors(IN SV* svn, IN SV* svk = 0)
PREINIT:
int status;
UV n, k, i, ndivisors, *divs;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
k = n;
if (status == 1 && svk != 0) {
status = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (k > n) k = n;
}
if (status != 1)
DISPATCHPP_RETURN();
if (GIMME_V == G_VOID) {
/* Nothing */
} else if (GIMME_V == G_SCALAR && k >= n) {
ndivisors = divisor_sum(n, 0);
PUSHs(sv_2mortal(newSVuv( ndivisors )));
} else {
divs = divisor_list(n, &ndivisors, k);
if (GIMME_V == G_SCALAR) {
PUSHs(sv_2mortal(newSVuv( ndivisors )));
} else {
EXTEND(SP, (EXTEND_TYPE)ndivisors);
for (i = 0; i < ndivisors; i++)
PUSHs(sv_2mortal(newSVuv( divs[i] )));
}
Safefree(divs);
}
void
trial_factor(IN SV* svn, ...)
ALIAS:
fermat_factor = 1
holf_factor = 2
squfof_factor = 3
lehman_factor = 4
prho_factor = 5
cheb_factor = 6
pplus1_factor = 7
pbrent_factor = 8
pminus1_factor = 9
ecm_factor = 10
PREINIT:
int nstatus, a1status, a2status, a3status;
UV n, arg1, arg2, arg3;
static const UV default_arg1[] =
{0, 64000000, 8000000, 4000000, 1, 4000000, 0, 200, 4000000, 1000000, 4000};
/* Trial, Fermat, Holf, SQUFOF, Lmn, PRHO, Cheb, P+1, Brent, P-1, ECM64 */
PPCODE:
if (ix == 10 && items > 2 && !SvOK(ST(1)) && SvOK(ST(2)))
croak("ecm_factor: B1 must be specified if B2 is specified");
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (items > 1) a1status=_validate_and_set(&arg1, aTHX_ ST(1), IFLAG_NONNEG);
else { a1status = 1; arg1 = default_arg1[ix]; }
if (items > 2) a2status=_validate_and_set(&arg2, aTHX_ ST(2), IFLAG_NONNEG);
else { a2status = 1; arg2 = 0; }
if (items > 3) a3status=_validate_and_set(&arg3, aTHX_ ST(3), IFLAG_NONNEG);
else { a3status = 1; arg3 = 0; }
/* See if we can do strint trial factoring here */
if (ix == 0 && !_XS_get_callgmp() && nstatus == 0 && a2status == 1
&& a1status == 1 && arg1 > 0 && arg1 <= 10000000) {
UV limit = (arg1 < 5) ? 5 : arg1;
STRLEN slen, cofactor_len;
const char *str = SvPV_nomg(svn, slen);
char *cofactor_str;
UV *fac_buf, cofactor_uv;
int nret = 0, nf, fi, nf_alloc = (int)(slen * 4) + 16;
New(0, fac_buf, nf_alloc, UV);
New(0, cofactor_str, slen + 1, char);
nf = strint_trial_factor(cofactor_str, &cofactor_len, &cofactor_uv,
fac_buf, str, slen, 2, (uint32_t)limit);
/* printf("strint trial through %u, found %d factors\n",
(uint32_t)limit, nf); */
EXTEND(SP, nf + 1);
for (fi = 0; fi < nf; fi++) {
PUSHs(sv_2mortal(newSVuv(fac_buf[fi])));
nret++;
}
if (nf == 0) {
PUSH_SV_CANONICAL(svn);
nret++;
} else if (cofactor_len == 0) {
if (cofactor_uv > 1) {
PUSHs(sv_2mortal(newSVuv(cofactor_uv)));
nret++;
}
} else {
PUSH_STR_CANONICAL(cofactor_str, cofactor_len);
nret++;
}
Safefree(fac_buf);
Safefree(cofactor_str);
XSRETURN(nret);
}
if (nstatus==0 || a1status==0 || a2status==0 || a3status==0 ||
(ix==10 && (_XS_get_callgmp() || !HAS_ECM64)))
DISPATCHPP_RETURN();
if (n == 0) XSRETURN_UV(0);
/* Small factors */
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 (is_prime(n)) { XPUSHs(sv_2mortal(newSVuv( n ))); }
else {
UV factors[MPU_MAX_FACTORS+1];
int i, nfactors = 0;
if (items < 3) {
if (ix == 8) arg2 = 1; /* use X*X + 1 */
if (ix == 9) arg2 = 10*arg1; /* B2 = 10*B1 */
}
if (ix == 10) arg3 = (items < 4) ? 100 : arg3; /* 100 curves */
switch (ix) {
case 0: nfactors = trial_factor (n, factors, 2, arg1); break;
case 1: nfactors = fermat_factor (n, factors, arg1); break;
case 2: nfactors = holf_factor (n, factors, arg1); break;
case 3: nfactors = squfof_factor (n, factors, arg1); break;
case 4: nfactors = lehman_factor (n, factors, arg1); break;
case 5: nfactors = prho_factor (n, factors, arg1); break;
case 6: nfactors = cheb_factor (n, factors, arg1, arg2); break;
case 7: nfactors = pplus1_factor (n, factors, arg1); break;
case 8: nfactors = pbrent_factor (n, factors, arg1, arg2); break;
case 9: nfactors = pminus1_factor (n, factors, arg1, arg2); break;
#if HAS_ECM64
case 10: nfactors = tinyecm64_factor(n, factors, arg1, arg2, arg3, 0); break;
#endif
default: break;
}
EXTEND(SP, (EXTEND_TYPE)nfactors);
for (i = 0; i < nfactors; i++)
PUSHs(sv_2mortal(newSVuv( factors[i] )));
}
void
divisor_sum(IN SV* svn, ...)
PREINIT:
UV n, k, sigma;
PPCODE:
if (items == 1) {
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
sigma = divisor_sum(n, 1);
if (n <= 1 || sigma != 0)
XSRETURN_UV(sigma);
}
} else {
SV* svk = ST(1);
if ( (!SvROK(svk) || (SvROK(svk) && SvTYPE(SvRV(svk)) != SVt_PVCV)) &&
_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG) ) {
sigma = divisor_sum(n, k);
if (n <= 1 || sigma != 0)
XSRETURN_UV(sigma);
}
}
DISPATCHPP_RETURN();
void aliquot_sum(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
UV sum = aliquot_sum(n);
if (n <= 1 || sum != 0)
XSRETURN_UV(sum);
}
DISPATCHPP_RETURN();
void abundance(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
UV sum = aliquot_sum(n);
if (n <= 1 || sum != 0)
XSRETURN_IV((IV)(sum-n));
}
DISPATCHPP_RETURN();
void sopf(IN SV* svn)
ALIAS:
sopfr = 1
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
UV sum = ix ? sopfr(n) : sopf(n);
XSRETURN_UV(sum);
}
DISPATCHPP_RETURN();
void prime_signature(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
factored_t nf = factorint(n);
uint32_t i, j, nfactors = nf.nfactors;
/* Insertion sort the exponents in descending order */
for (i = 1; i < nfactors; i++) {
uint8_t t = nf.e[i];
for (j = i; j > 0 && nf.e[j-1] < t; j--)
nf.e[j] = nf.e[j-1];
nf.e[j] = t;
}
if (GIMME_V == G_ARRAY) {
EXTEND(SP, (EXTEND_TYPE)nfactors);
for (i = 0; i < nfactors; i++)
PUSH_NPARITY(nf.e[i]);
XSRETURN(nfactors);
} else {
UV S = nfactors > 0, p = 0;
for (i = 0; i < nfactors; i++)
{ p = next_prime(p); S *= ipow(p,nf.e[i]); }
XSRETURN_UV(S);
}
}
DISPATCHPP_RETURN();
void jordan_totient(IN SV* svk, IN SV* svn)
PREINIT:
int kstatus, nstatus;
UV k, n;
PPCODE:
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (kstatus != 1)
croak("jordan_totient: k must fit in a UV");
if (nstatus == 1) {
UV ret = jordan_totient(k, n);
if (ret != UV_MAX)
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void powersum(IN SV* svn, IN SV* svk)
PREINIT:
int kstatus, nstatus;
UV k, n;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (kstatus != 1)
croak("powersum: k must fit in a UV");
if (nstatus == 1) {
UV ret = powersum(n, k);
if (ret != UV_MAX)
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void
ramanujan_sum(IN SV* sva, IN SV* svn)
ALIAS:
legendre_phi = 1
smooth_count = 2
rough_count = 3
PREINIT:
int astatus, nstatus;
UV a, n, ret;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_NONNEG);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (astatus != 0 && nstatus != 0) {
switch (ix) {
case 0: if (a < 1 || n < 1) XSRETURN_IV(0);
{
UV g = a / gcd_ui(a,n);
int m = moebius(g);
if (m == 0 || a == g) RETURN_NPARITY(m);
XSRETURN_IV( m * (totient(a) / totient(g)) );
}
break;
case 1: ret = legendre_phi(a, n); break;
case 2: ret = debruijn_psi(a, n); break;
case 3:
default: ret = buchstab_phi(a, n); break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void almost_prime_count(IN SV* svk, IN SV* svn)
ALIAS:
almost_prime_count_approx = 1
almost_prime_count_lower = 2
almost_prime_count_upper = 3
omega_prime_count = 4
PREINIT:
UV k, n, ret;
PPCODE:
if (_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG) &&
_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
k < BITS_PER_WORD) {
ret = 0;
switch (ix) {
case 0: ret = almost_prime_count(k, n); break;
case 1: ret = almost_prime_count_approx(k, n); break;
case 2: ret = almost_prime_count_lower(k, n); break;
case 3: ret = almost_prime_count_upper(k, n); break;
case 4: ret = omega_prime_count(k, n); break;
default: break;
}
XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void nth_almost_prime(IN SV* svk, IN SV* svn)
ALIAS:
nth_almost_prime_approx = 1
nth_almost_prime_lower = 2
nth_almost_prime_upper = 3
PREINIT:
UV k, n, max;
PPCODE:
if (_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG) &&
_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
k < BITS_PER_WORD) {
UV ret = 0;
if (n == 0 || (k == 0 && n > 1)) XSRETURN_UNDEF;
max = max_almost_prime_count(k);
if (max > 0 && n <= max) {
switch (ix) {
case 0: ret = nth_almost_prime(k, n); break;
case 1: ret = nth_almost_prime_approx(k, n); break;
case 2: ret = nth_almost_prime_lower(k, n); break;
case 3: ret = nth_almost_prime_upper(k, n); break;
}
if (ret != 0) XSRETURN_UV(ret);
}
}
DISPATCHPP_RETURN();
void nth_omega_prime(IN SV* svk, IN SV* svn)
PREINIT:
UV k, n, max, ret;
PPCODE:
if (_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG) &&
_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) &&
k < 16) {
if (n == 0 || (k == 0 && n > 1)) XSRETURN_UNDEF;
max = max_omega_prime_count(k);
if (max > 0 && n <= max) {
ret = nth_omega_prime(k, n);
if (ret != 0) XSRETURN_UV(ret);
}
}
DISPATCHPP_RETURN();
void powmod(IN SV* sva, IN SV* svg, IN SV* svn)
ALIAS:
rootmod = 1
PREINIT:
int astatus, gstatus, nstatus, retundef;
UV a, g, n, ret;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
gstatus = _validate_and_set(&g, aTHX_ svg, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && gstatus != 0 && nstatus != 0) {
if (n == 0) XSRETURN_UNDEF;
if (n == 1) XSRETURN_UV(0);
_mod_with(&a, astatus, n);
retundef = ret = 0;
if (ix == 0) {
retundef = !prep_pow_inv(&a,&g,gstatus,n);
if (!retundef) ret = powmod(a, g, n);
} else {
retundef = !(prep_pow_inv(&a,&g,gstatus,n) && rootmod(&ret,a,g,n));
}
if (retundef) XSRETURN_UNDEF;
XSRETURN_UV(ret);
}
if (!_XS_get_callgmp() && ix == 0) {
STRLEN lena, leng, lenn, rlen;
const char *sa = SvPV_nomg(sva, lena), *sg = SvPV_nomg(svg, leng), *sn = SvPV_nomg(svn, lenn);
SV* tmp = sv_2mortal(newSV(0 + lenn));
rlen = strint_powmod(SvPVX(tmp), sa,lena, sg,leng, sn,lenn);
if (rlen > 0)
RETURN_PVSV_CANONICAL(tmp, rlen);
}
DISPATCHPP_RETURN();
void addmod(IN SV* sva, IN SV* svb, IN SV* svn)
ALIAS:
submod = 1
mulmod = 2
divmod = 3
znlog = 4
PREINIT:
int astatus, bstatus, nstatus, retundef;
UV a, b, n, ret;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
bstatus = _validate_and_set(&b, aTHX_ svb, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && bstatus != 0 && nstatus != 0) {
if (n == 0) XSRETURN_UNDEF;
if (n == 1) XSRETURN_UV(0);
_mod_with(&a, astatus, n);
_mod_with(&b, bstatus, n);
retundef = ret = 0;
switch (ix) {
case 0: ret = addmod(a, b, n); break;
case 1: ret = submod(a, b, n); break;
case 2: ret = mulmod(a, b, n); break;
case 3: b = modinverse(b, n);
if (b == 0) retundef = 1;
else ret = mulmod(a, b, n);
break;
case 4: if (b == 0) { /* g == 0: sequence is 1, 0, 0, 0, ... */
if (a == 1) { ret = 0; }
else if (a == 0) { ret = 1; }
else { retundef = 1; }
} else {
ret = znlog(a, b, n);
if (ret == 0 && a != 1) retundef = 1;
}
break;
default: break;
}
if (retundef) XSRETURN_UNDEF;
XSRETURN_UV(ret);
}
#if HAVE_FACTOR128
if (ix <= 2) {
uint128_t a128, b128, n128, ret128;
int asign, bsign;
if (xs_sv_to_uint128_signmag(aTHX_ &a128, &asign, sva) &&
xs_sv_to_uint128_signmag(aTHX_ &b128, &bsign, svb) &&
xs_sv_to_uint128_abs(aTHX_ &n128, svn)) {
if (n128 == 0) XSRETURN_UNDEF;
if (n128 == 1) XSRETURN_UV(0);
a128 = mod_with128(a128, asign, n128);
b128 = mod_with128(b128, bsign, n128);
switch (ix) {
case 0: ret128 = muladdmod128_s(1, a128, b128, n128, 0); break;
case 1: ret128 = muladdmod128_s(1, a128, b128, n128, 1); break;
default: ret128 = muladdmod128_s(a128, b128, 0, n128, 0); break;
}
RETURN_U128(ret128);
}
}
#endif
if (!_XS_get_callgmp() && ix <= 2) {
STRLEN lena, lenb, lenn, rlen;
const char *sa = SvPV_nomg(sva, lena), *sb = SvPV_nomg(svb, lenb), *sn = SvPV_nomg(svn, lenn);
SV* tmp = sv_2mortal(newSV(0 + lenn));
switch (ix) {
case 0: rlen = strint_muladdmod_s(SvPVX(tmp), "1",1, sa,lena, sb,lenb, 0, sn,lenn); break;
case 1: rlen = strint_muladdmod_s(SvPVX(tmp), "1",1, sa,lena, sb,lenb, 1, sn,lenn); break;
default: rlen = strint_muladdmod_s(SvPVX(tmp), sa,lena, sb,lenb, "0",1, 0, sn,lenn); break;
}
RETURN_PVSV_CANONICAL(tmp, rlen);
}
DISPATCHPP_RETURN();
void muladdmod(IN SV* sva, IN SV* svb, IN SV* svc, IN SV* svn)
ALIAS:
mulsubmod = 1
PREINIT:
int astatus, bstatus, cstatus, nstatus;
UV a, b, c, n, ret;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
bstatus = _validate_and_set(&b, aTHX_ svb, IFLAG_ANY);
cstatus = _validate_and_set(&c, aTHX_ svc, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && bstatus != 0 && cstatus != 0 && nstatus != 0) {
if (n == 0) XSRETURN_UNDEF;
if (n == 1) XSRETURN_UV(0);
_mod_with(&a, astatus, n);
_mod_with(&b, bstatus, n);
_mod_with(&c, cstatus, n);
ret = (ix==0) ? muladdmod(a,b,c,n) : mulsubmod(a,b,c,n);
XSRETURN_UV(ret);
}
#if HAVE_FACTOR128
{
uint128_t a128, b128, c128, n128, ret128;
int asign, bsign, csign;
if (xs_sv_to_uint128_signmag(aTHX_ &a128, &asign, sva) &&
xs_sv_to_uint128_signmag(aTHX_ &b128, &bsign, svb) &&
xs_sv_to_uint128_signmag(aTHX_ &c128, &csign, svc) &&
xs_sv_to_uint128_abs(aTHX_ &n128, svn)) {
if (n128 == 0) XSRETURN_UNDEF;
if (n128 == 1) XSRETURN_UV(0);
a128 = mod_with128(a128, asign, n128);
b128 = mod_with128(b128, bsign, n128);
c128 = mod_with128(c128, csign, n128);
ret128 = muladdmod128_s(a128, b128, c128, n128, ix);
RETURN_U128(ret128);
}
}
#endif
if (!_XS_get_callgmp()) {
STRLEN lena, lenb, lenc, lenn, rlen;
const char *sa = SvPV_nomg(sva, lena), *sb = SvPV_nomg(svb, lenb), *sc = SvPV_nomg(svc, lenc), *sn = SvPV_nomg(svn, lenn);
SV* tmp = sv_2mortal(newSV(0 + lenn));
rlen = strint_muladdmod_s(SvPVX(tmp), sa,lena, sb,lenb, sc,lenc, ix, sn,lenn);
RETURN_PVSV_CANONICAL(tmp,rlen);
}
DISPATCHPP_RETURN();
void binomialmod(IN SV* svn, IN SV* svk, IN SV* svm)
PREINIT:
int nstatus, kstatus, mstatus;
UV ret, n, k, m;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_ANY);
mstatus = _validate_and_set(&m, aTHX_ svm, IFLAG_ABS);
if (nstatus != 0 && kstatus != 0 && mstatus != 0) {
if (m == 0) XSRETURN_UNDEF;
if (m == 1) XSRETURN_UV(0);
if ( (nstatus == 1 && (kstatus == -1 || k > n)) ||
(nstatus ==-1 && (kstatus == -1 && k > n)) )
XSRETURN_UV(0);
if (kstatus == -1) k = n - k;
if (nstatus == -1) {
UV nabs = neg_iv(n);
if (k == 0) XSRETURN_UV(1);
if (nabs > UV_MAX - k + 1)
DISPATCHPP_RETURN();
n = nabs + k - 1;
}
if (binomialmod(&ret, n, k, m)) {
if ((nstatus == -1) && (k & 1) && ret != 0) ret = m-ret;
XSRETURN_UV(ret);
}
}
DISPATCHPP_RETURN();
void factorialmod(IN SV* sva, IN SV* svn)
PREINIT:
int astatus, nstatus;
UV a, n;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_NONNEG);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && nstatus != 0) {
if (n == 0) XSRETURN_UNDEF;
if (n == 1) XSRETURN_UV(0);
XSRETURN_UV( factorialmod(a, n) );
}
DISPATCHPP_RETURN();
void invmod(IN SV* sva, IN SV* svn)
ALIAS:
znorder = 1
sqrtmod = 2
negmod = 3
PREINIT:
int astatus, nstatus;
UV a, n, r, retok;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && nstatus != 0) {
if (n == 0) XSRETURN_UNDEF;
if (n == 1) XSRETURN_UV((ix==1) ? 1 : 0); /* znorder different */
_mod_with(&a, astatus, n);
retok = r = 0;
switch (ix) {
case 0: retok = r = modinverse(a, n); break;
case 1: retok = r = znorder(a, n); break;
case 2: retok = sqrtmod(&r, a, n); break;
case 3:
default: retok = 1; r = negmod(a, n); break;
}
if (retok == 0) XSRETURN_UNDEF;
XSRETURN_UV(r);
}
DISPATCHPP_RETURN();
void allsqrtmod(IN SV* sva, IN SV* svn)
PREINIT:
int astatus, nstatus;
UV a, n, i, numr, *roots;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && nstatus != 0) {
if (n == 0) RETURN_NOTHING();
_mod_with(&a, astatus, n);
roots = allsqrtmod(&numr, a, n);
if (roots != 0) {
if (GIMME_V != G_ARRAY) {
PUSHs(sv_2mortal(newSVuv(numr)));
} else {
EXTEND(SP, (EXTEND_TYPE)numr);
for (i = 0; i < numr; i++)
PUSHs(sv_2mortal(newSVuv(roots[i])));
}
Safefree(roots);
} else if (GIMME_V != G_ARRAY) {
PUSHs(sv_2mortal(newSVuv(0)));
}
} else {
DISPATCHPP_RETURN();
}
void allrootmod(IN SV* sva, IN SV* svg, IN SV* svn)
PREINIT:
int astatus, gstatus, nstatus;
UV a, g, n, i, numr, *roots;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
gstatus = _validate_and_set(&g, aTHX_ svg, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && gstatus != 0 && nstatus != 0) {
if (n == 0) RETURN_NOTHING();
if (n == 1) XSRETURN_UV(GIMME_V == G_ARRAY ? 0 : 1);
_mod_with(&a, astatus, n);
if (!prep_pow_inv(&a,&g,gstatus,n))
RETURN_NOTHING();
roots = allrootmod(&numr, a, g, n);
if (roots != 0) {
if (GIMME_V != G_ARRAY) {
PUSHs(sv_2mortal(newSVuv(numr)));
} else {
EXTEND(SP, (EXTEND_TYPE)numr);
for (i = 0; i < numr; i++)
PUSHs(sv_2mortal(newSVuv(roots[i])));
}
Safefree(roots);
} else if (GIMME_V != G_ARRAY) {
PUSHs(sv_2mortal(newSVuv(0)));
}
} else {
DISPATCHPP_RETURN();
}
void is_primitive_root(IN SV* sva, IN SV* svn)
PREINIT:
int astatus, nstatus;
UV a, n;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && nstatus != 0) {
if (n == 0) XSRETURN_UNDEF;
_mod_with(&a, astatus, n);
RETURN_NPARITY( is_primitive_root(a,n,0) );
}
DISPATCHPP_RETURN();
void qnr(IN SV* svn)
ALIAS:
znprimroot = 1
PREINIT:
UV n, r;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_ABS)) {
if (n == 0) XSRETURN_UNDEF;
if (ix == 0) {
r = qnr(n);
} else {
r = znprimroot(n);
if (r == 0 && n != 1) XSRETURN_UNDEF;
}
if (r < 100) RETURN_NPARITY(r);
else XSRETURN_UV(r);
}
DISPATCHPP_RETURN();
void
is_smooth(IN SV* svn, IN SV* svk)
ALIAS:
is_rough = 1
PREINIT:
UV n, k;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_ABS) &&
_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG)) {
RETURN_NPARITY( (ix == 0) ? is_smooth(n,k) : is_rough(n,k) );
}
DISPATCHPP_RETURN();
void
is_omega_prime(IN SV* svk, IN SV* svn)
ALIAS:
is_almost_prime = 1
PREINIT:
UV n, k;
int nstatus, kstatus;
PPCODE:
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (kstatus != 0 && nstatus != 0) {
int res = (nstatus != 1) ? 0
: (ix == 0) ? is_omega_prime(k, n)
: is_almost_prime(k, n);
RETURN_NPARITY(res);
}
#if HAVE_FACTOR128
if (kstatus != 0 && nstatus == 0) {
factored128_t nf;
if (xs_factorintp128_sv(aTHX_ &nf, svn)) {
uint32_t nfac = (ix == 0) ? factored128p_distinct_factors(&nf)
: factored128p_total_factors(&nf);
RETURN_NPARITY(nfac == k);
}
}
#endif
DISPATCHPP_RETURN();
void is_divisible(IN SV* svn, IN SV* svd, ...)
PREINIT:
UV n, d, ret;
size_t i;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_ABS) &&
_validate_and_set(&d, aTHX_ svd, IFLAG_ABS)) {
int status = 1;
ret = d==0 ? (n==0) : n % d == 0;
for (i = 2; i < (size_t)items && !ret; i++) {
if ((status = _validate_and_set(&d, aTHX_ ST(i), IFLAG_ABS)) != 1)
break;
ret = d==0 ? (n==0) : n % d == 0;
}
if (status == 1) RETURN_NPARITY(ret);
}
DISPATCHPP_RETURN();
void is_congruent(IN SV* svn, IN SV* svc, IN SV* svd)
PREINIT:
UV n, c, d;
int nstatus, cstatus, dstatus;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
cstatus = _validate_and_set(&c, aTHX_ svc, IFLAG_ANY);
dstatus = _validate_and_set(&d, aTHX_ svd, IFLAG_ABS);
if (nstatus != 0 && cstatus != 0 && dstatus != 0) {
if (d != 0) {
_mod_with(&n, nstatus, d);
_mod_with(&c, cstatus, d);
}
RETURN_NPARITY( n == c );
}
DISPATCHPP_RETURN();
void valuation(IN SV* svn, IN SV* svk)
PREINIT:
UV n, k;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_ABS) &&
_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG)) {
if (k <= 1) croak("valuation: k must be > 1");
if (n == 0) XSRETURN_UNDEF;
RETURN_NPARITY(valuation(n, k));
}
DISPATCHPP_RETURN();
void remove_factors(IN SV* svn, IN SV* svk)
ALIAS:
remove_factors_exp = 1
PREINIT:
UV n, k, e;
int nstatus, kstatus;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (nstatus != 0 && kstatus != 0) {
if (k <= 1) croak("%s: k must be > 1", SUBNAME);
if (n == 0) {
XPUSHs(&PL_sv_undef);
if (ix == 1) XPUSHs(&PL_sv_undef);
} else {
if (nstatus == -1) n = (UV)neg_iv(n);
e = valuation_remainder(n, k, &n);
XPUSHs(sv_2mortal( nstatus == 1 ? newSVuv(n) : newSViv(neg_iv(n)) ));
if (ix == 1) XPUSHs(sv_2mortal(newSVuv(e)));
}
} else {
DISPATCHPP_RETURN();
}
void is_powerful(IN SV* svn, IN SV* svk = 0);
ALIAS:
powerful_count = 1
sumpowerful = 2
nth_powerful = 3
PREINIT:
int nstatus;
UV n, ret, k = 2;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, (ix < 3) ? IFLAG_ANY: IFLAG_NONNEG);
if (nstatus != 0 && (!svk || _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG))) {
if (nstatus == -1) RETURN_NPARITY(0);
if (ix == 0) RETURN_NPARITY( is_powerful(n, k) );
if (ix == 1) XSRETURN_UV( powerful_count(n, k) );
if (ix == 2) {
if (n == 0) XSRETURN_UV(0);
ret = sumpowerful(n, k);
} else {
if (n == 0) XSRETURN_UNDEF;
ret = nth_powerful(n, k);
}
/* ret=0: nth_powerful / sumpowerful result > UV_MAX, so go to PP/GMP */
if (ret > 0) XSRETURN_UV(ret);
}
DISPATCHPP_RETURN();
void kronecker(IN SV* sva, IN SV* svb)
PREINIT:
int k;
PPCODE:
if (xs_kronecker_result(aTHX_ sva, svb, &k))
RETURN_NPARITY( k );
DISPATCHPP_RETURN();
void is_qr(IN SV* sva, IN SV* svn)
PREINIT:
int astatus, nstatus;
UV a, n;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (astatus != 0 && nstatus != 0) {
if (n == 0) XSRETURN_UNDEF;
if (n == 1) RETURN_NPARITY(1);
_mod_with(&a, astatus, n);
RETURN_NPARITY( is_qr(a,n) );
}
DISPATCHPP_RETURN();
void addint(IN SV* sva, IN SV* svb)
ALIAS:
subint = 1
mulint = 2
divint = 3
modint = 4
cdivint = 5
powint = 7
PREINIT:
SV* slowret;
int astatus, bstatus, is_uv;
UV uvret;
IV ivret;
PPCODE:
if (addint_try_native_result(aTHX_ ix, sva, svb, SUBNAME,
&astatus, &bstatus, &is_uv, &uvret, &ivret)) {
if (is_uv) XSRETURN_UV(uvret);
XSRETURN_IV(ivret);
}
slowret = addint_try_slow_result(aTHX_ ix, sva, svb, astatus, bstatus, SUBNAME);
if (slowret != NULL)
RETURN_SV(slowret);
/* Others get dispatched here */
DISPATCHPP_RETURN();
void muladdint(IN SV* sva, IN SV* svb, IN SV* svc)
ALIAS:
mulsubint = 1
PREINIT:
int astatus, bstatus, cstatus;
UV a, b, c;
PPCODE:
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_ANY);
bstatus = _validate_and_set(&b, aTHX_ svb, IFLAG_ANY);
cstatus = _validate_and_set(&c, aTHX_ svc, IFLAG_ANY);
if (astatus != 0 && bstatus != 0 && cstatus != 0) {
if (astatus == -1) a = neg_iv(a);
if (bstatus == -1) b = neg_iv(b);
if (cstatus == -1) c = neg_iv(c);
{
int sign;
UV hi, lo;
if (muladd_uv_signmag(&sign, &hi, &lo, a, b, c,
astatus, bstatus, (ix == 1) ? -cstatus : cstatus))
RETURN_SIGNMAG_UV_UV(sign, hi, lo);
}
}
{
STRLEN lena, lenb, lenc, rlen;
const char *sa = SvPV_nomg(sva, lena), *sb = SvPV_nomg(svb, lenb), *sc = SvPV_nomg(svc, lenc);
SV* tmp = sv_2mortal(newSV(2 + lena + lenb + lenc));
rlen = strint_muladd_s(SvPVX(tmp), sa,lena, sb,lenb, sc,lenc, ix);
RETURN_PVSV_CANONICAL(tmp,rlen);
}
DISPATCHPP_RETURN();
void add1int(IN SV* svn)
ALIAS:
sub1int = 1
PREINIT:
SV* slowret;
int status, is_uv;
UV uvret;
IV ivret;
PPCODE:
if (add1_try_native_result(aTHX_ ix, svn, &status, &is_uv, &uvret, &ivret)) {
if (is_uv) XSRETURN_UV(uvret);
XSRETURN_IV(ivret);
}
slowret = add1_try_slow_result(aTHX_ ix, svn);
if (slowret != NULL)
RETURN_SV(slowret);
DISPATCHPP_RETURN();
void absint(IN SV* svn)
ALIAS:
negint = 1
PREINIT:
UV n;
PPCODE:
if (ix == 0) {
if (_validate_and_set(&n, aTHX_ svn, IFLAG_ABS))
XSRETURN_UV(n);
} else {
int status = _validate_and_set(&n, aTHX_ svn, IFLAG_IV);
if (status == -1) XSRETURN_UV(neg_iv(n));
else if (status == 1) XSRETURN_IV(neg_iv(n));
}
TRY_MAGIC_UNARY(svn, ix==0 ? abs_amg : neg_amg);
{
STRLEN len;
const char* s = SvPV_nomg(svn, len);
SV* tmp = sv_2mortal(newSV(1 + len));
len = ix==0 ? strint_abs(SvPVX(tmp),s,len) : strint_neg(SvPVX(tmp),s,len);
RETURN_PVSV_CANONICAL(tmp,len);
}
DISPATCHPP_RETURN();
void signint(IN SV* svn)
ALIAS:
is_odd = 1
is_even = 2
PPCODE:
RETURN_NPARITY(xs_sign_parity_result(aTHX_ svn, SUBNAME, ix));
void cmpint(IN SV* sva, IN SV* svb)
PREINIT:
int ret;
PPCODE:
ret = xs_cmpint_result(aTHX_ sva, svb);
RETURN_NPARITY(ret);
void logint(IN SV* svn, IN SV* svb, IN SV* svret = 0)
PREINIT:
UV n, b;
int nstatus, bstatus;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_POS);
bstatus = _validate_and_set(&b, aTHX_ svb, IFLAG_POS);
if (bstatus != 0 && b < 2) croak("logint: base must be > 1");
if (items >= 3 && !xs_is_sv_scalar_ref(svret))
croak("logint: third argument not a scalar reference");
if (nstatus != 0 && bstatus != 0) {
UV ilog = logint(n, b);
if (svret) sv_setuv(SvRV(svret), ipow(b,ilog));
XSRETURN_UV(ilog);
}
/* GMP backend is fast. PP not so fast even with Math::GMPz. */
/* Fast path: strint_logint for bigint n and UV base */
if (bstatus == 1 && _XS_get_callgmp() < 47) {
STRLEN lenn;
const char* sn = SvPV_nomg(svn, lenn);
UV result = strint_logint(sn, lenn, b);
if (result != UV_MAX) {
if (svret) {
STRLEN blen;
const char* bstr = SvPV_nomg(sv_2mortal(newSVuv(b)), blen);
SV* vtmp = sv_2mortal(newSV(lenn + 2));
STRLEN vlen = strint_pow(SvPVX(vtmp), bstr, blen, result, lenn+2);
SvCUR_set(vtmp, vlen); SvPOK_on(vtmp); *SvEND(vtmp) = '\0';
sv_setsv(SvRV(svret), xs_to_canonical(aTHX_ vtmp));
}
XSRETURN_UV(result);
}
}
DISPATCHPP_RETURN_GMPIF(svret == 0 && (bstatus || _XS_get_callgmp() >= 54));
void rootint(IN SV* svn, IN SV* svk, IN SV* svret = 0)
PREINIT:
UV n, k;
int nstatus, kstatus;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_POS);
if (items >= 3 && !xs_is_sv_scalar_ref(svret))
croak("rootint: third argument not a scalar reference");
if (nstatus != 0 && kstatus != 0) {
UV iroot = rootint(n, k);
if (svret) sv_setuv(SvRV(svret), ipow(iroot,k));
XSRETURN_UV(iroot);
}
/* Fast path: strint_rootint for bigint n, UV k */
if (kstatus == 1 && _XS_get_callgmp() < 40) {
STRLEN lenn;
const char* sn = SvPV_nomg(svn, lenn);
SV* tmp = sv_2mortal(newSV(lenn + 2));
STRLEN rlen = strint_rootint(SvPVX(tmp), sn, lenn, k);
if (rlen > 0) {
if (svret) {
SV* vtmp = sv_2mortal(newSV(lenn + 2));
STRLEN vlen = strint_pow(SvPVX(vtmp), SvPVX(tmp), rlen, k, lenn + 2);
SvCUR_set(vtmp, vlen); SvPOK_on(vtmp); *SvEND(vtmp) = '\0';
sv_setsv(SvRV(svret), xs_to_canonical(aTHX_ vtmp));
}
RETURN_PVSV_CANONICAL(tmp, rlen);
}
}
DISPATCHPP_RETURN_GMPIF(svret == 0 && (kstatus || _XS_get_callgmp() >= 54));
void crootint(IN SV* svn, IN SV* svk)
PREINIT:
UV n, k;
int nstatus, kstatus;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_POS);
if (nstatus != 0 && kstatus != 0)
XSRETURN_UV(crootint(n, k));
DISPATCHPP_RETURN();
void divrem(IN SV* sva, IN SV* svb)
ALIAS:
fdivrem = 1
cdivrem = 2
tdivrem = 3
PREINIT:
int astatus, bstatus;
UV D, d;
IV iD, id;
PPCODE:
astatus = _validate_and_set(&D, aTHX_ sva, IFLAG_ANY);
bstatus = _validate_and_set(&d, aTHX_ svb, IFLAG_ANY);
if (astatus != 0 && bstatus != 0 && d == 0)
croak("%s: divide by zero", SUBNAME);
if (astatus == 1 && bstatus == 1 && (ix != 2 || D % d == 0)) {
XPUSHs(sv_2mortal(newSVuv( D / d )));
XPUSHs(sv_2mortal(newSVuv( D % d )));
XSRETURN(2);
} else if (ix == 2 && astatus == 1 && bstatus == 1 && d <= (UV)IV_MAX) {
/* Exact division was handled above */
XPUSHs(sv_2mortal(newSVuv( D/d + 1 )));
XPUSHs(sv_2mortal(newSViv( ((IV)D%d) - d )));
XSRETURN(2);
} else if (astatus != 0 && bstatus != 0 &&
_validate_and_set((UV*)&iD, aTHX_ sva, IFLAG_IV) != 0 &&
_validate_and_set((UV*)&id, aTHX_ svb, IFLAG_IV) != 0 &&
iD != IV_MIN && id != IV_MIN) {
/* Both values, and their negations, fit in an IV */
IV q, r;
switch (ix) {
case 0: edivrem(&q, &r, iD, id); break;
case 1: fdivrem(&q, &r, iD, id); break;
case 2: cdivrem(&q, &r, iD, id); break;
case 3:
default: tdivrem(&q, &r, iD, id); break;
}
XPUSHs(sv_2mortal(newSViv( q )));
XPUSHs(sv_2mortal(newSViv( r )));
XSRETURN(2);
}
DISPATCHPP_RETURN();
void lshiftint(IN SV* svn, IN SV* svk = 0)
ALIAS:
rshiftint = 1
rashiftint = 2
PREINIT:
int nstatus, kstatus, nix;
UV n, k, nk;
PPCODE:
nix = ix;
if (items == 1) {
kstatus = 1;
k = 1;
} else {
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_ANY);
if (kstatus == -1) {
k = neg_iv(k);
nix = !ix; /* 0 => 1, 1 => 0, 2 => 0 */
}
}
if (kstatus != 0) {
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (k == 0)
XSRETURN(1);
if (nstatus != 0 && nix > 0 && k >= BITS_PER_WORD) /* Big right shift */
XSRETURN_IV(nstatus == -1 && nix==2 ? -1 : 0);
if (nstatus == 1 && k < BITS_PER_WORD) {
if (nix > 0) XSRETURN_UV(n >> k); /* Right shift */
if ( ((n << k) >> k) == n) XSRETURN_UV(n << k); /* Left shift */
/* Fall through -- left shift needs more bits */
} else if (nstatus == -1 && nix > 0 && k < BITS_PER_WORD) {
n = neg_iv(n);
nk = n >> k;
XSRETURN_IV( nix == 1 ? -nk : (nk<<k)==n ? -nk : -nk-1 );
} else if (nstatus == -1 && nix == 0 && k+1 < BITS_PER_WORD) {
n = neg_iv(n);
nk = n << k;
if ((nk << 1) >> (k+1) == n)
XSRETURN_IV(-nk);
/* Fall through -- left shift needs more bits */
}
if (k < 1600) { /* String fast path for relatively small shifts */
STRLEN len;
const char* s = SvPV_nomg(svn, len);
SV* tmp = sv_2mortal(newSV(len + (nix == 0 ? 3 + k*.302 : 0)));
switch (nix) {
case 0: len = strint_lshiftint( SvPVX(tmp), s, len, k); break;
case 1: len = strint_rshiftint( SvPVX(tmp), s, len, k); break;
default: len = strint_rashiftint(SvPVX(tmp), s, len, k); break;
}
if (len > 0) RETURN_PVSV_CANONICAL(tmp,len);
}
}
DISPATCHPP_RETURN();
void
gcdext(IN SV* sva, IN SV* svb)
PREINIT:
IV u, v, d, a, b;
PPCODE:
if (_validate_and_set((UV*)&a, aTHX_ sva, IFLAG_IV) &&
_validate_and_set((UV*)&b, aTHX_ svb, IFLAG_IV) &&
a != IV_MIN && b != IV_MIN) {
d = gcdext(a, b, &u, &v, 0, 0);
XPUSHs(sv_2mortal(newSViv( u )));
XPUSHs(sv_2mortal(newSViv( v )));
XPUSHs(sv_2mortal(newSViv( d )));
} else {
DISPATCHPP_RETURN();
}
void
stirling(IN SV* svn, IN SV* svm, IN UV type = 1)
PREINIT:
UV n, m;
int nstatus, mstatus;
PPCODE:
if (type != 1 && type != 2 && type != 3)
croak("stirling: type must be 1, 2, or 3");
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
mstatus = _validate_and_set(&m, aTHX_ svm, IFLAG_NONNEG);
if (nstatus == 1 && mstatus == 1) {
if (n == m)
XSRETURN_UV(1);
if (n == 0 || m == 0 || m > n)
XSRETURN_UV(0);
if (type == 3) {
UV s = stirling3(n, m);
if (s != 0) XSRETURN_UV(s);
} else if (type == 2) {
IV s = stirling2(n, m);
if (s != 0) XSRETURN_IV(s);
} else if (type == 1) {
IV s = stirling1(n, m);
if (s != 0) XSRETURN_IV(s);
}
}
DISPATCHPP_RETURN();
NV
_XS_ExponentialIntegral(IN SV* x)
ALIAS:
_XS_LogarithmicIntegral = 1
_XS_RiemannZeta = 2
_XS_RiemannR = 3
_XS_LambertW = 4
PREINIT:
NV nv, ret;
CODE:
nv = !SvROK(x) ? SvNV(x) : STRTONV(SvPV_nolen(x));
switch (ix) {
case 0: ret = Ei(nv); break;
case 1: ret = Li(nv); break;
case 2: ret = (NV) ld_riemann_zeta(nv); break;
case 3: ret = (NV) RiemannR(nv,0); break;
case 4:
default:ret = lambertw(nv); break;
}
RETVAL = ret;
OUTPUT:
RETVAL
void euler_phi(IN SV* svlo, IN SV* svhi = 0)
ALIAS:
moebius = 1
PREINIT:
UV lo, hi;
int lostatus, histatus;
PPCODE:
lostatus = _validate_and_set(&lo, aTHX_ svlo, IFLAG_ANY);
if (items == 1) {
if (lostatus != 0 && ix == 0)
XSRETURN_UV(lostatus == -1 ? 0 : totient(lo));
if (lostatus != 0 && ix == 1)
RETURN_NPARITY(moebius(lostatus == -1 ? neg_iv(lo) : lo));
#if HAVE_FACTOR128
if (ix == 1) {
uint128_t n128;
if (xs_sv_to_uint128_abs(aTHX_ &n128, svlo))
RETURN_NPARITY(moebius128(n128));
}
#endif
DISPATCHPP_RETURN(); /* single argument dispatch */
}
histatus = _validate_and_set(&hi, aTHX_ svhi, IFLAG_ANY);
if (GIMME_V != G_ARRAY && lostatus == 1 && histatus == 1 && hi != UV_MAX)
XSRETURN_UV(lo > hi ? 0 : hi - lo + 1);
if (GIMME_V != G_ARRAY) { /* Scalar context = return range */
STRLEN lolen, hilen, rlen;
const char *lostr = SvPV_nomg(svlo, lolen);
const char *histr = SvPV_nomg(svhi, hilen);
SV *tmp;
if (strint_cmp(lostr, lolen, histr, hilen) > 0)
XSRETURN_UV(0);
tmp = sv_2mortal(newSV((lolen > hilen ? lolen : hilen) + 2));
rlen = strint_sub(SvPVX(tmp), histr, hilen, lostr, lolen); /* hi-lo */
rlen = strint_add(SvPVX(tmp), SvPVX(tmp), rlen, "1", 1); /* +1 */
RETURN_PVSV_CANONICAL(tmp, rlen);
}
if (lostatus != 1 || histatus != 1)
DISPATCHPP_RETURN(); /* ranged list return */
if (lo > hi) XSRETURN_EMPTY;
{
UV i, count;
bool appendmax = (hi == UV_MAX); /* Strip hi=UV_MAX from the loop */
if (appendmax) hi--;
count = hi-lo+1;
if (count > MAX_EXTEND - appendmax)
croak("%s: range too large for list return", SUBNAME);
EXTEND(SP, (EXTEND_TYPE)count + appendmax);
if (count > 0) {
if (ix == 0) {
UV arrlo = (lo < 100) ? 0 : lo;
UV *totients = range_totient(arrlo, hi);
for (i = 0; i < count; i++)
PUSHs(sv_2mortal(newSVuv(totients[i+lo-arrlo])));
Safefree(totients);
} else {
signed char* mu;
UV seglo, seghi;
void* mctx = start_segment_moebius(lo, hi, &mu);
while (next_segment_moebius(mctx, &seglo, &seghi)) {
for (i = seglo; i <= seghi; i++)
PUSH_NPARITY(mu[i-seglo]);
}
end_segment_moebius(mctx);
}
}
if (appendmax) {
if (ix == 0) PUSHs(sv_2mortal(newSVuv(totient(UV_MAX))));
else PUSH_NPARITY(-1); /* moebius 2^{32,64,128}-1 = -1 */
}
}
void
dedekind_psi(IN SV* svn)
PREINIT:
int nstatus;
UV n, r;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (nstatus != 0) {
if (nstatus == -1) XSRETURN_UV(0);
r = dedekind_psi(n);
if (n == 0 || r > 0) XSRETURN_UV(r);
}
DISPATCHPP_RETURN();
void sqrtint(IN SV* svn)
ALIAS:
carmichael_lambda = 1
exp_mangoldt = 2
PREINIT:
UV n, r;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
r = 0;
switch (ix) {
case 0: r = isqrt(n); break;
case 1: r = carmichael_lambda(n); break;
case 2: r = exp_mangoldt(n); break;
default: break;
}
XSRETURN_UV(r);
}
/* GMP backend is very fast. PP sorta-fast with Math::GMPz/Math::GMP */
/* Fast path: strint_rootint for bigint n */
if (ix == 0 && _XS_get_callgmp() < 40) {
STRLEN lenn;
const char* sn = SvPV_nomg(svn, lenn);
SV* tmp = sv_2mortal(newSV(lenn + 2));
STRLEN rlen = strint_rootint(SvPVX(tmp), sn, lenn, 2);
if (rlen > 0) RETURN_PVSV_CANONICAL(tmp,rlen);
}
DISPATCHPP_RETURN();
void hammingweight(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_ABS))
XSRETURN_UV(popcnt(n));
if (_XS_get_callgmp() < 47) {
char* ptr; STRLEN len; ptr = SvPV(svn, len);
XSRETURN_UV(mpu_popcount_string(ptr, len));
}
DISPATCHPP_RETURN();
void prime_omega(IN SV* svn)
ALIAS:
prime_bigomega = 1
is_square_free = 2
PREINIT:
UV n, ret;
int status;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (status != 0) {
ret = 0;
switch (ix) {
case 0: ret = prime_omega(n); break;
case 1: ret = prime_bigomega(n); break;
case 2: ret = is_square_free(n); break;
default: break;
}
RETURN_NPARITY(ret);
}
#if HAVE_FACTOR128
{
uint128_t n128;
if (xs_sv_to_uint128_abs(aTHX_ &n128, svn)) {
if (ix < 2) {
factored128_t nf;
factorintp128(&nf, n128);
ret = (ix == 0) ? factored128p_distinct_factors(&nf)
: factored128p_total_factors(&nf);
} else ret = (moebius128(n128) != 0);
RETURN_NPARITY(ret);
}
}
#endif
DISPATCHPP_RETURN();
void pisano_period(IN SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
UV r = pisano_period(n);
if (r != UV_MAX)
XSRETURN_UV(r);
}
DISPATCHPP_RETURN();
void factorial(IN SV* svn)
ALIAS:
subfactorial = 1
bell_number = 2
fubini = 3
primorial = 4
pn_primorial = 5
catalan_number = 6
partitions = 7
partitionsq = 8
consecutive_integer_lcm = 9
PREINIT:
UV n, r = 0;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG) != 1)
croak("%s: n must fit in native unsigned integer", SUBNAME);
switch(ix) {
case 0: r = factorial(n); break;
case 1: r = subfactorial(n); break;
case 2: r = bell_number(n); break;
case 3: r = fubini(n); break;
case 4: r = primorial(n); break;
case 5: r = pn_primorial(n); break;
case 6: r = catalan_number(n); break;
case 7: r = npartitions(n); break;
case 8: r = npartitionsq(n); break;
case 9: r = consecutive_integer_lcm(n); break;
default: break;
}
if (n == 0 || r > 0)
XSRETURN_UV(r);
DISPATCHPP_RETURN();
void integer_complexity(IN SV* svn)
PREINIT:
int nstatus;
UV n, r;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (nstatus == 0 || n > (UV)IV_MAX)
croak("integer_complexity: n must fit in native signed integer");
r = integer_complexity(n); /* Make sure to call it even if n=0 */
if (n == 0) XSRETURN_UNDEF;
XSRETURN_UV(r);
void sumtotient(IN SV* svn)
PREINIT:
UV n, r;
#if HAVE_FACTOR128 && HAVE_SUMTOTIENT128
uint64_t n64;
uint128_t sum128;
#endif
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
r = sumtotient(n);
if (n == 0 || r > 0) XSRETURN_UV(r);
#if HAVE_FACTOR128 && HAVE_SUMTOTIENT128
if (sumtotient128(n, &sum128)) /* 64-bit overflowed, try 128-bit. */
RETURN_U128(sum128);
#endif
}
#if HAVE_FACTOR128 && HAVE_SUMTOTIENT128
else if (xs_sv_to_uint64(aTHX_ &n64, svn)) {
if (sumtotient128(n64, &sum128))
RETURN_U128(sum128);
}
#endif
DISPATCHPP_RETURN();
void binomial(IN SV* svn, IN SV* svk)
PREINIT:
int nstatus, kstatus;
UV n, k, ret;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_ANY);
if (nstatus != 0 && kstatus != 0) {
if ( (nstatus == 1 && (kstatus == -1 || k > n)) ||
(nstatus ==-1 && (kstatus == -1 && k > n)) )
XSRETURN_UV(0);
if (kstatus == -1) {
k = n - k; /* n<0,k<=n: (-1)^(n-k) * binomial(-k-1,n-k) */
kstatus = 1;
}
if (nstatus == -1) {
UV nabs = neg_iv(n);
if (k == 0) XSRETURN_UV(1);
if (nabs <= UV_MAX - k + 1) {
UV ntop = nabs + k - 1;
ret = binomial(ntop, k);
if (ret > 0 && ret <= (UV)IV_MAX)
XSRETURN_IV( (IV)ret * ((k&1) ? -1 : 1) );
/* The result overflowed. Use strint. */
if (ntop <= UINT32_MAX && k <= UINT32_MAX && _XS_get_callgmp() < 53) {
STRLEN rlen;
char *rstr = strint_binomial_u32_u32((uint32_t)ntop, (uint32_t)k, &rlen);
if (rstr)
RETURN_SIGN_STRINT_STR((k&1) ? -1 : 1, rstr, rlen);
}
}
} else if (nstatus == 1) {
ret = binomial(n, k);
if (ret != 0) XSRETURN_UV(ret);
/* The result overflowed. Use strint. */
if (n <= UINT32_MAX && k <= UINT32_MAX && _XS_get_callgmp() < 53) {
STRLEN rlen;
char *rstr = strint_binomial_u32_u32((uint32_t)n, (uint32_t)k, &rlen);
if (rstr)
RETURN_SIGN_STRINT_STR(1, rstr, rlen);
}
}
}
if (kstatus == 1 && k <= UINT32_MAX && _XS_get_callgmp() < 53) {
STRLEN snlen, rlen;
const char *sn;
char *rstr;
SvGETMAGIC(svn);
sn = SvPV_nomg(svn, snlen);
rstr = strint_binomial_u32(sn, snlen, (uint32_t)k, &rlen);
if (rstr)
RETURN_SIGN_STRINT_STR(1, rstr, rlen);
}
DISPATCHPP_RETURN();
void multifactorial(IN SV* svn, IN SV* svk)
PREINIT:
int nstatus, kstatus;
UV n, k, r;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_POS);
if (nstatus == 1 && kstatus == 1) {
r = multifactorial(n, k);
if (n == 0 || r > 0) XSRETURN_UV(r);
}
DISPATCHPP_RETURN();
void falling_factorial(IN SV* svn, IN SV* svk)
ALIAS:
rising_factorial = 1
PREINIT:
int nstatus, kstatus;
UV n, k;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_IV);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);
if (nstatus == 1 && kstatus == 1) {
UV ret = (ix==0) ? falling_factorial(n,k) : rising_factorial(n,k);
if (ret != UV_MAX) XSRETURN_UV(ret);
} else if (nstatus == -1 && kstatus == 1) {
IV in = (IV)n;
IV ret = (ix==0) ? falling_factorial_s(in,k) : rising_factorial_s(in,k);
if (ret != IV_MAX) XSRETURN_IV(ret);
}
DISPATCHPP_RETURN();
void floor_sum(IN SV* svn, IN SV* svm, IN SV* sva, IN SV* svb)
PREINIT:
int nstatus, mstatus, astatus, bstatus;
UV n, m, a, b;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
mstatus = _validate_and_set(&m, aTHX_ svm, IFLAG_POS);
astatus = _validate_and_set(&a, aTHX_ sva, IFLAG_NONNEG);
bstatus = _validate_and_set(&b, aTHX_ svb, IFLAG_NONNEG);
if (nstatus == 1 && mstatus == 1 && astatus == 1 && bstatus == 1) {
UV sum = floor_sum(n,m,a,b);
if (sum != UV_MAX)
XSRETURN_UV(sum);
}
DISPATCHPP_RETURN();
void mertens(IN SV* svlo, IN SV* svhi = 0)
PREINIT:
UV lo = 1, hi;
PPCODE:
if ((items == 1 && _validate_and_set(&hi, aTHX_ svlo, IFLAG_NONNEG)) ||
(items == 2 && _validate_and_set(&lo, aTHX_ svlo, IFLAG_NONNEG) &&
_validate_and_set(&hi, aTHX_ svhi, IFLAG_NONNEG))) {
RETURN_NPARITY(mertens_range(lo, hi));
}
DISPATCHPP_RETURN();
void liouville(IN SV* svn)
ALIAS:
sumliouville = 1
is_pillai = 2
is_congruent_number = 3
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
IV r = 0;
switch(ix) {
case 0: r = liouville(n); break;
case 1: r = sumliouville(n); break;
case 2: r = pillai_v(n); break;
case 3: r = is_congruent_number(n); break;
default: break;
}
RETURN_NPARITY(r);
}
#if HAVE_FACTOR128
if (ix == 0) {
uint128_t n128;
if (xs_sv_to_uint128(aTHX_ &n128, svn)) {
factored128_t nf;
factorintp128(&nf, n128);
RETURN_NPARITY((factored128p_total_factors(&nf) & 1) ? -1 : 1);
}
}
#endif
DISPATCHPP_RETURN();
void hclassno(IN SV* svn)
PREINIT:
UV n, r;
int status;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == -1)
XSRETURN_IV(0);
if (status == 1) {
if (n == 0)
XSRETURN_IV(-1);
r = hclassno(n);
if (r != UV_MAX)
XSRETURN_UV(r);
}
DISPATCHPP_RETURN();
void ramanujan_tau(IN SV* svn)
PREINIT:
UV n;
int status;
PPCODE:
status = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (status == -1)
XSRETURN_IV(0);
if (status == 1) {
IV r = ramanujan_tau(n);
if (n == 0 || r != 0)
RETURN_NPARITY(r);
}
DISPATCHPP_RETURN();
int _is_congruent_number_filter(IN UV n)
CODE:
RETVAL = is_congruent_number_filter(n);
OUTPUT:
RETVAL
bool _is_congruent_number_tunnell(IN UV n)
CODE:
RETVAL = is_congruent_number_tunnell(n);
OUTPUT:
RETVAL
void chebyshev_theta(IN SV* svn)
ALIAS:
chebyshev_psi = 1
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
NV r = (ix==0) ? chebyshev_theta(n) : chebyshev_psi(n);
XSRETURN_NV(r);
}
DISPATCHPP_RETURN(); /* Result is FP */
#define RETURN_SET_REF(ps) /* Return sorted set values (ps is iset_t pointer) */ \
{ \
UV *sdata; \
size_t slen = iset_size(ps); \
int sign = iset_sign(ps); \
New(0, sdata, slen, UV); \
iset_allvals(ps, sdata); \
iset_destroy(ps); \
RETURN_LIST_REF( slen, sdata, sign ); \
}
#define RETURN_EMPTY_SET_REF() RETURN_EMPTY_LIST_REF()
void sumset(IN SV* sva, IN SV* svb = 0)
PROTOTYPE: $;$
PREINIT:
int atype, btype, stype, sign;
UV *ra, *rb;
size_t alen, blen, i, j;
iset_t s;
PPCODE:
atype = arrayref_to_int_array(aTHX_ &alen, &ra, 1, sva, "sumset arg 1");
if (svb == 0 || atype == IARR_TYPE_BAD) {
rb = ra;
blen = alen;
btype = atype;
} else {
btype = arrayref_to_int_array(aTHX_ &blen, &rb, 1, svb, "sumset arg 2");
}
if (alen == 0 || blen == 0) {
if (rb != ra) Safefree(rb);
Safefree(ra);
RETURN_EMPTY_SET_REF();
}
if (atype == IARR_TYPE_BAD || btype == IARR_TYPE_BAD)
stype = IARR_TYPE_BAD;
else
stype = type_of_sumset(atype, btype, ra[0],ra[alen-1], rb[0],rb[blen-1]);
if (stype == IARR_TYPE_BAD) {
if (rb != ra) Safefree(rb);
Safefree(ra);
DISPATCHPP_RETURN();
}
sign = IARR_TYPE_TO_STATUS(stype);
/* Sumset */
iset_create(&s, 10UL * (alen+blen));
for (i = 0; i < alen; i++)
for (j = 0; j < blen; j++)
iset_add(&s, ra[i]+rb[j], sign);
if (rb != ra) Safefree(rb);
Safefree(ra);
RETURN_SET_REF(&s);
void setbinop(IN SV* block, IN SV* sva, IN SV* svb = 0)
PROTOTYPE: &$;$
PREINIT:
int atype, btype;
UV *ra, *rb;
Size_t alen, blen;
CODE:
/* Must be CODE and not PPCODE */
atype = arrayref_to_int_array(aTHX_ &alen, &ra, 1, sva, "setbinop arg 1");
if (svb == 0 || atype == IARR_TYPE_BAD) {
rb = ra;
blen = alen;
btype = atype;
} else {
btype = arrayref_to_int_array(aTHX_ &blen, &rb, 1, svb, "setbinop arg 2");
}
if (alen == 0 || blen == 0) {
if (rb != ra) Safefree(rb);
Safefree(ra);
RETURN_EMPTY_SET_REF();
}
if (atype != IARR_TYPE_BAD && btype != IARR_TYPE_BAD) {
iset_t s;
Size_t i, j;
GV *agv, *bgv;
SV *asv, *bsv;
UV ret;
CV *subcv;
int status = 0;
SETSUBREF(subcv, block);
agv = gv_fetchpv("a", GV_ADD, SVt_PV);
bgv = gv_fetchpv("b", GV_ADD, SVt_PV);
asv = NEWSVINT(0,0);
bsv = NEWSVINT(0,0);
SAVEFREESV(asv);
SAVEFREESV(bsv);
SAVESPTR(GvSV(agv));
SAVESPTR(GvSV(bgv));
GvSV(agv) = asv;
GvSV(bgv) = bsv;
iset_create(&s, 4UL * ((size_t)alen + (size_t)blen + 2));
/* Native multicall pre-5.10.1 cannot safely redispatch with these args.*/
#if USE_MULTICALL && (!defined(REAL_MULTICALL) || PERL_VERSION_GE(5,10,1))
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
SV *cbret;
I32 gimme = G_SCALAR;
SC_PUSH_MULTICALL(subcv);
for (i = 0; i < alen; i++) {
for (j = 0; j < blen; j++) {
FASTSETSVINT(asv, atype == IARR_TYPE_POS, ra[i]);
FASTSETSVINT(bsv, btype == IARR_TYPE_POS, rb[j]);
SC_MULTICALL_SCALAR(cbret);
status = _validate_and_set(&ret, aTHX_ cbret, IFLAG_ANY);
if (status != 0) iset_add(&s, ret, status);
if (status == 0 || iset_is_invalid(&s)) break;
}
if (j < blen) break;
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
{
for (i = 0; i < alen; i++) {
for (j = 0; j < blen; j++) {
SV* cbret;
FASTSETSVINT(asv, atype == IARR_TYPE_POS, ra[i]);
FASTSETSVINT(bsv, btype == IARR_TYPE_POS, rb[j]);
cbret = xs_call_cv_noinput_1_sv(aTHX_ subcv);
status = _validate_and_set(&ret, aTHX_ cbret, IFLAG_ANY);
if (status != 0) iset_add(&s, ret, status);
if (status == 0 || iset_is_invalid(&s)) break;
}
if (j < blen) break;
}
}
if (status != 0 && !iset_is_invalid(&s)) {
if (rb != ra) Safefree(rb);
Safefree(ra);
RETURN_SET_REF(&s);
}
iset_destroy(&s);
}
if (rb != ra) Safefree(rb);
Safefree(ra);
DISPATCHPP_RETURN();
void setunion(IN SV* sva, IN SV* svb)
PROTOTYPE: $$
ALIAS:
setintersect = 1
setminus = 2
setdelta = 3
PREINIT:
SV *ret;
PPCODE:
/* Alias ix values intentionally match set_op_t. */
if (xs_set_op(aTHX_ sva, svb, (set_op_t)ix, &ret, SUBNAME)) {
ST(0) = sv_2mortal(ret);
XSRETURN(1);
}
DISPATCHPP_RETURN();
void set_is_disjoint(IN SV* sva, IN SV* svb)
PROTOTYPE: $$
ALIAS:
set_is_equal = 1
set_is_subset = 2
set_is_proper_subset = 3
set_is_superset = 4
set_is_proper_superset = 5
set_is_proper_intersection = 6
PREINIT:
int ret;
PPCODE:
/* Alias ix values intentionally match set_relation_op_t. */
if (xs_set_relation(aTHX_ sva, svb, (set_relation_op_t)ix, &ret, SUBNAME))
RETURN_NPARITY(ret);
DISPATCHPP_RETURN();
void setcontains(IN SV* sva, ...)
ALIAS:
setcontainsany = 1
PROTOTYPE: $@
PREINIT:
UV b;
AV *ava;
int bstatus, subset, findall;
Size_t alen, blen, i;
DECL_ARREF(arb);
PPCODE:
CHECK_ARRAYREF(sva); /* First argument is a set as array ref */
ava = (AV*) SvRV(sva);
alen = av_count(ava);
if (items < 2)
RETURN_NPARITY(ix == 0 ? 1 : 0);
if (SvMAGICAL(ava) || !AvREAL(ava)) /* Punt these to Perl */
DISPATCHPP_RETURN();
findall = ix == 0 ? 1 : 0;
if (items == 2 && SvROK(ST(1)) && SvTYPE(SvRV(ST(1))) == SVt_PVAV) {
set_cache_t svcache;
USE_ARREF(arb, ST(1), SUBNAME, AR_READ);
/* If setcontainsany and B is bigger than A, swap them for performance. */
if (ix == 1 && len_arb > alen && svarr_arb != 0) {
ava = avp_arb;
alen = len_arb;
USE_ARREF(arb, ST(0), SUBNAME, AR_READ);
}
blen = len_arb;
subset = ix == 0 && blen > alen ? 0 : findall;
_sc_init_cache(&svcache);
/* setcontains: if we find anything that is NOT in SETA, return 0
* setcontainsany: if we find anything that IS in SETA, return 1 */
for (i = 0; i < blen && subset == findall; i++) {
bstatus = _validate_and_set(&b, aTHX_ FETCH_ARREF(arb,i), IFLAG_ANY);
subset = is_in_set(aTHX_ ava, &svcache, bstatus, b);
REFRESH_ARREF(arb);
}
} else {
UV *rb;
int btype = array_to_int_array(aTHX_ &blen, &rb, 1, &ST(1), items-1);
bstatus = IARR_TYPE_TO_STATUS(btype);
subset = bstatus == 0 ? -1 : ix == 0 && blen > alen ? 0 : findall;
if (blen <= 4) {
for (i = 0; i < blen && subset == findall; i++)
subset = is_in_set(aTHX_ ava, 0, bstatus, rb[i]);
} else {
set_cache_t svcache;
_sc_init_cache(&svcache);
for (i = 0; i < blen && subset == findall; i++)
subset = is_in_set(aTHX_ ava, &svcache, bstatus, rb[i]);
}
Safefree(rb);
}
if (subset != -1)
RETURN_NPARITY(subset);
DISPATCHPP_RETURN();
void setinsert(IN SV* sva, ...)
PROTOTYPE: $@
PREINIT:
AV *ava;
Size_t alen, blen, i, j;
UV *rb;
int btype, bstatus;
PPCODE:
CHECK_ARRAYREF(sva); /* First argument is a set as array ref */
ava = (AV*) SvRV(sva);
alen = av_count(ava);
if (items < 2)
RETURN_NPARITY(0);
CHECK_AV_NOT_READONLY(ava); /* We intend to modify it */
if (SvMAGICAL(ava) || !AvREAL(ava)) /* Punt these to Perl */
DISPATCHPP_RETURN();
if (SvROK(ST(1)) && SvTYPE(SvRV(ST(1))) == SVt_PVAV) {
if (items != 2)
croak("setinsert: expected integer list or single array reference");
btype = arrayref_to_int_array(aTHX_ &blen, &rb, 1, ST(1), "setinsert");
} else {
btype = array_to_int_array(aTHX_ &blen, &rb, 1, &ST(1), items-1);
}
bstatus = IARR_TYPE_TO_STATUS(btype);
if (bstatus != 0 && blen <= 4) {
int res = 0;
size_t nins = 0;
for (i = 0; res >= 0 && i < blen; i++) {
res = ins_into_set(aTHX_ ava, bstatus, rb[i]);
nins += (res > 0);
}
if (res >= 0) { Safefree(rb); RETURN_NPARITY(nins); }
} else if (bstatus != 0) {
size_t nbeg, nmid, nend, nmidcheck;
int alostatus, ahistatus;
UV alo, ahi;
set_cache_t svcache;
/* 1. ava is empty. push everything and we're done. */
if (alen == 0) {
av_extend(ava, blen);
for (i = 0; i < blen; i++)
av_push(ava, NEWSVINT(bstatus, rb[i]));
Safefree(rb);
RETURN_NPARITY(blen);
}
_sc_init_cache(&svcache);
/* Get hi and lo values of set. */
if (_sc_set_bounds(aTHX_ ava, &svcache, alen,
&alostatus, &ahistatus, &alo, &ahi) >= 0) {
if (_sign_cmp(alostatus,alo,ahistatus,ahi) > 0) {
Safefree(rb);
croak("%s: expected numerically ascending sorted input", SUBNAME);
}
/* Both lo/hi are not bigint, so there are no bigints in the set. */
nbeg = nend = nmid = 0;
/* 1. Find out how many elements go in front. */
while (nbeg < blen && _sign_cmp(bstatus,rb[nbeg],alostatus,alo) < 0)
nbeg++;
/* 2. Find out how many elements go at the end. */
while (nend < blen-nbeg && _sign_cmp(bstatus,rb[blen-1-nend],ahistatus,ahi) > 0)
nend++;
/* 3. In-place insert everything in the middle. */
nmidcheck = blen - nbeg - nend;
if (nmidcheck > 0) {
size_t *insert_idx;
croak("%s: expected sorted input, found bigint value in interior", SUBNAME);
}
if (index > 0) {
insert_sv[nmid] = NEWSVINT(bstatus,rb[i]); /* Value to insert */
insert_idx[nmid] = (Size_t)index-1; /* Where to insert */
nmid++;
}
}
av_extend(ava, alen + nmid + nbeg + nend);
if (nmid > 0) {
SV** arr;
Size_t index_lastorig = alen-1;
Size_t index_moveto = index_lastorig + nmid;
/* Push new values on end so Perl calculates array correctly. */
for (i = 0; i < nmid; i++)
av_push(ava, insert_sv[i]);
arr = AvARRAY(ava);
/* SV* pointer manipulation to insert new values in place. */
for (i = 0; i < nmid; i++) {
size_t j = nmid-1-i;
size_t idx = insert_idx[j];
size_t nmove = index_lastorig - idx + 1;
if (nmove > 0) {
size_t moveto = index_moveto - nmove + 1;
memmove(arr+moveto, arr+idx, sizeof(SV*) * nmove);
index_lastorig -= nmove;
index_moveto -= nmove;
}
arr[index_moveto--] = insert_sv[j];
}
}
Safefree(insert_sv);
Safefree(insert_idx);
}
/* 4. Insert at front */
if (nbeg > 0) {
av_unshift(ava, nbeg);
for (i = 0; i < nbeg; i++)
av_store(ava, i, NEWSVINT(bstatus, rb[i]));
}
/* 5. Push onto back */
if (nend > 0) {
for (i = 0; i < nend; i++)
av_push(ava, NEWSVINT(bstatus, rb[blen-nend+i]));
}
Safefree(rb);
RETURN_NPARITY(nbeg+nmid+nend);
}
}
Safefree(rb);
DISPATCHPP_RETURN();
void setremove(IN SV* sva, ...)
PROTOTYPE: $@
PREINIT:
AV *ava;
Size_t alen, blen, i;
UV *rb;
int btype, bstatus;
PPCODE:
CHECK_ARRAYREF(sva); /* First argument is a set as array ref */
ava = (AV*) SvRV(sva);
alen = av_count(ava);
if (alen == 0 || items < 2)
RETURN_NPARITY(0);
CHECK_AV_NOT_READONLY(ava); /* We intend to modify it */
if (SvMAGICAL(ava) || !AvREAL(ava)) /* Punt these to Perl */
DISPATCHPP_RETURN();
if (SvROK(ST(1)) && SvTYPE(SvRV(ST(1))) == SVt_PVAV) {
if (items != 2)
croak("setremove: expected integer list or single array reference");
btype = arrayref_to_int_array(aTHX_ &blen, &rb, 1, ST(1), "setremove");
} else {
btype = array_to_int_array(aTHX_ &blen, &rb, 1, &ST(1), items-1);
}
if (btype != IARR_TYPE_BAD) {
bstatus = IARR_TYPE_TO_STATUS(btype);
if (blen <= 5 || alen <= 20) { /* SIMPLE DELETE LOOP */
int res = 0;
size_t ndel = 0;
for (i = 0; res >= 0 && i < blen; i++) {
res = del_from_set(aTHX_ ava, bstatus, rb[i]);
if (res > 0) ndel++;
}
if (res >= 0) { Safefree(rb); RETURN_NPARITY(ndel); }
} else if (blen < 500 || (blen*100) < alen) { /* ONE PASS DELETE */
Size_t *del_idx, ndel = 0;
set_cache_t svcache;
_sc_init_cache(&svcache);
/* Create index list to remove */
New(0, del_idx, blen, Size_t);
for (i = 0; i < blen; i++) {
SSize_t index = index_in_set(aTHX_ ava, &svcache, bstatus, rb[i]);
if (index < 0) {
Safefree(del_idx);
Safefree(rb);
croak("%s: expected sorted input, found bigint value in interior", SUBNAME);
}
if (index > 0)
del_idx[ndel++] = (Size_t)index-1;
}
Safefree(rb);
if (ndel > 0) {
SV **arr = AvARRAY(ava);
size_t to = del_idx[0];
for (i = 0; i < ndel; i++) {
size_t idx = del_idx[i];
size_t beg = idx+1;
size_t len = (i+1) >= ndel ? alen-beg : del_idx[i+1]-beg;
SvREFCNT_dec_NN(arr[idx]);
if (len > 0) {
memmove(arr+to, arr+beg, sizeof(SV*) * len);
to += len;
}
}
Zero(arr + alen - ndel, ndel, SV*);
av_fill(ava, alen-ndel-1);
}
Safefree(del_idx);
RETURN_NPARITY(ndel);
} else { /* CLEAR AND GREP */
int atype, astatus, del_complete = 0;
UV *ra = 0;
atype = arrayref_to_int_array(aTHX_ &alen, &ra, 1, sva, SUBNAME);
if (CAN_COMBINE_IARR_TYPES(atype,btype)) {
size_t ia = 0, ib = 0;
int pcmp = (atype == IARR_TYPE_NEG || btype == IARR_TYPE_NEG) ? 0 : 1;
astatus = IARR_TYPE_TO_STATUS(atype);
av_clear(ava);
while (ia < alen && ib < blen) {
if (ra[ia] == rb[ib]) {
ia++; ib++;
} else {
if (SIGNED_CMP_LT(pcmp, ra[ia], rb[ib])) av_push(ava, NEWSVINT(astatus, ra[ia++]));
else ib++;
}
}
while (ia < alen) av_push(ava, NEWSVINT(astatus, ra[ia++]));
del_complete = 1;
}
Safefree(ra);
Safefree(rb);
if (del_complete) RETURN_NPARITY(alen - av_count(ava));
}
}
DISPATCHPP_RETURN();
void setinvert(IN SV* sva, ...)
PROTOTYPE: $@
PREINIT:
AV *ava;
Size_t alen, blen, i;
UV *rb;
int btype, bstatus;
PPCODE:
CHECK_ARRAYREF(sva);
ava = (AV*) SvRV(sva);
alen = av_count(ava);
if (items < 2)
RETURN_NPARITY(0);
CHECK_AV_NOT_READONLY(ava);
if (SvMAGICAL(ava) || !AvREAL(ava))
DISPATCHPP_RETURN();
if (SvROK(ST(1)) && SvTYPE(SvRV(ST(1))) == SVt_PVAV) {
if (items != 2)
croak("setinvert: expected integer list or single array reference");
btype = arrayref_to_int_array(aTHX_ &blen, &rb, 1, ST(1), "setinvert");
} else {
btype = array_to_int_array(aTHX_ &blen, &rb, 1, &ST(1), items-1);
}
if (btype != IARR_TYPE_BAD) {
if (blen == 0) {
Safefree(rb);
RETURN_NPARITY(0);
}
bstatus = IARR_TYPE_TO_STATUS(btype);
if (blen <= 4 || alen <= 20) { /* SIMPLE TOGGLE LOOP */
IV ndelta = 0;
int res = 0;
for (i = 0; res >= 0 && i < blen; i++) {
res = del_from_set(aTHX_ ava, bstatus, rb[i]);
if (res > 0) { ndelta--; } /* found and removed */
else if (res == 0) { /* not found, insert */
res = ins_into_set(aTHX_ ava, bstatus, rb[i]);
if (res > 0) ndelta++;
}
}
if (res >= 0) {
Safefree(rb);
ST(0) = sv_2mortal(newSViv(ndelta));
XSRETURN(1);
}
} else { /* MERGE-STYLE SYMMETRIC DIFFERENCE */
int atype, astatus, done = 0;
UV *ra = 0;
Size_t old_alen = alen;
atype = arrayref_to_int_array(aTHX_ &alen, &ra, 1, sva, SUBNAME);
if (CAN_COMBINE_IARR_TYPES(atype, btype)) {
size_t ia = 0, ib = 0;
int pcmp = (atype == IARR_TYPE_NEG || btype == IARR_TYPE_NEG) ? 0 : 1;
astatus = IARR_TYPE_TO_STATUS(atype);
av_clear(ava);
while (ia < alen && ib < blen) {
if (ra[ia] == rb[ib]) { ia++; ib++; }
else if (SIGNED_CMP_LT(pcmp, ra[ia], rb[ib])) av_push(ava, NEWSVINT(astatus, ra[ia++]));
else av_push(ava, NEWSVINT(bstatus, rb[ib++]));
}
while (ia < alen) av_push(ava, NEWSVINT(astatus, ra[ia++]));
while (ib < blen) av_push(ava, NEWSVINT(bstatus, rb[ib++]));
done = 1;
}
Safefree(ra);
if (done) {
Safefree(rb);
ST(0) = sv_2mortal(newSViv((IV)av_count(ava) - (IV)old_alen));
XSRETURN(1);
}
}
}
Safefree(rb);
DISPATCHPP_RETURN();
void is_sidon_set(IN SV* sva)
PROTOTYPE: $
PREINIT:
int ret;
PPCODE:
if (xs_is_sidon_set(aTHX_ sva, &ret))
RETURN_NPARITY(ret);
DISPATCHPP_RETURN();
void is_sumfree_set(IN SV* sva)
PROTOTYPE: $
PREINIT:
int ret;
PPCODE:
if (xs_is_sumfree_set(aTHX_ sva, &ret))
RETURN_NPARITY(ret);
DISPATCHPP_RETURN();
void toset(...)
PROTOTYPE: @
PREINIT:
int type;
size_t len;
UV *L;
PPCODE:
if (items == 0) RETURN_EMPTY_SET_REF();
if (items == 1 && SvROK(ST(0)) && SvTYPE(SvRV(ST(0))) == SVt_PVAV)
croak("toset: expected integer list, not array reference");
type = array_to_int_array(aTHX_ &len, &L, 1, &ST(0), items);
if (type != IARR_TYPE_BAD)
RETURN_LIST_REF(len, L, type != IARR_TYPE_NEG);
Safefree(L);
DISPATCHPP_RETURN();
void vecsort(...)
PROTOTYPE: @
ALIAS:
vecrsort = 1
PREINIT:
int type;
size_t len;
UV *L;
PPCODE:
if (items == 0)
RETURN_NOTHING();
if (SvROK(ST(0)) && SvTYPE(SvRV(ST(0))) == SVt_PVAV) {
if (items != 1)
croak("%s: expected integer list or single array reference", SUBNAME);
type = arrayref_to_int_array(aTHX_ &len, &L, 0, ST(0), SUBNAME);
} else {
type = array_to_int_array(aTHX_ &len, &L, 0, &ST(0), items);
}
if (GIMME_V != G_ARRAY) /* In scalar context, return number of elements */
XSRETURN_UV(len);
if (type == IARR_TYPE_ANY || type == IARR_TYPE_POS) {
sort_uv_array(L, len);
} else if (type == IARR_TYPE_NEG) {
sort_iv_array((IV*)L, len);
} else {
Safefree(L);
DISPATCHPP_RETURN();
}
if (ix)
reverse_uv_array(L, len);
RETURN_LIST_VALS( len, L, (type != IARR_TYPE_NEG) );
void vecsorti(IN SV* sva)
PROTOTYPE: $
ALIAS:
vecrsorti = 1
PREINIT:
int type;
size_t i, len;
UV *L;
SV **arr;
AV *ava;
PPCODE:
CHECK_ARRAYREF(sva);
ava = (AV*) SvRV(sva);
CHECK_AV_NOT_READONLY(ava); /* We intend to modify it */
if (SvMAGICAL(ava) || !AvREAL(ava)) /* Punt these to Perl */
DISPATCHPP_RETURN();
type = arrayref_to_int_array(aTHX_ &len, &L, 0, sva, SUBNAME);
/* If we really wanted to optimize small values, the reading function
* could create a mask like:
* mask |= (istatus == 1) ? n : (n ^ (n<<1));
* then we know if the input is 8-bit, 16-bit, 32-bit, etc.
*/
if (type == IARR_TYPE_ANY || type == IARR_TYPE_POS) {
sort_uv_array(L, len);
} else if (type == IARR_TYPE_NEG) {
sort_iv_array((IV*)L, len);
} else {
Safefree(L);
DISPATCHPP_RETURN();
}
arr = AvARRAY(ava);
for (i = 0; i < len; i++)
FASTSETSVINT(arr[i], type == IARR_TYPE_POS, L[ix ? len-i-1 : i]);
Safefree(L);
XSRETURN(1);
void numtoperm(IN SV* svn, IN SV* svk)
PREINIT:
UV k, n, fn;
int nstatus, kstatus;
int i, S[32];
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_ANY);
if (nstatus != 0 && kstatus != 0 && n < 32) {
if (n == 0)
RETURN_NOTHING();
fn = factorial(n);
if (fn != 0) {
_mod_with(&k, kstatus, fn);
if (num_to_perm(k, n, S)) {
if (GIMME_V != G_ARRAY) XSRETURN_UV(n);
EXTEND(SP, (EXTEND_TYPE)n);
for (i = 0; i < (int)n; i++)
PUSH_NPARITY( S[i] );
XSRETURN(n);
}
}
}
DISPATCHPP_RETURN();
void permtonum(IN SV* svp)
PREINIT:
UV val, num;
Size_t i, plen;
DECL_ARREF(avp);
PPCODE:
USE_ARREF(avp, svp, SUBNAME, AR_READ);
plen = len_avp;
if (plen <= 20) {
int V[21], A[21] = {0};
for (i = 0; i < plen; i++) {
SV *iv = FETCH_ARREF(avp,i);
int status = _validate_and_set(&val, aTHX_ iv, IFLAG_NONNEG);
REFRESH_ARREF(avp);
if (status != 1)
break;
if (val >= plen || A[val] != 0) break;
A[val] = i+1;
V[i] = val;
}
if (i >= plen && perm_to_num(plen, V, &num))
XSRETURN_UV(num);
}
DISPATCHPP_RETURN();
void randperm(IN SV* svn, IN SV* svk = 0)
PREINIT:
int nstatus, kstatus;
UV n, k, i, *S;
dMY_CXT;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (items == 1) { kstatus = nstatus; k = nstatus ? n : 0; }
else { kstatus = _validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG);}
if (nstatus == 0)
DISPATCHPP_RETURN();
if (kstatus == 0 || k > n)
k = n;
if (k > (UV)IV_MAX)
croak("randperm: k must fit in native signed integer");
if (GIMME_V != G_ARRAY) XSRETURN_IV(k);
if (k == 0) XSRETURN_EMPTY;
if (k > (UV)MAX_SSIZET || k > (UV)(MAX_SIZET / sizeof(UV)))
croak("randperm: requested permutation is too large");
New(0, S, k, UV);
randperm(MY_CXT.randcxt, n, k, S);
EXTEND(SP, (EXTEND_TYPE)k);
for (i = 0; i < k; i++) {
if (n < 2*CINTS) PUSH_NPARITY(S[i]);
else PUSHs(sv_2mortal(newSVuv(S[i])));
}
Safefree(S);
void shuffle(...)
PROTOTYPE: @
PREINIT:
SSize_t i, j;
void* randcxt;
dMY_CXT;
PPCODE:
if (GIMME_V != G_ARRAY) XSRETURN_IV(items);
if (items == 0) XSRETURN_EMPTY;
for (i = 0, randcxt = MY_CXT.randcxt; i < items-1; i++) {
j = urandomm64(randcxt, items-i);
{ SV* t = ST(i); ST(i) = ST(i+j); ST(i+j) = t; }
}
XSRETURN(items);
void vecsample(IN SV* svk, ...)
PROTOTYPE: $@
PREINIT:
void *randcxt;
UV k;
Size_t nitems, i;
dMY_CXT;
PPCODE:
if (_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG) != 1)
DISPATCHPP_RETURN();
if (items == 1 || k == 0)
RETURN_NOTHING();
randcxt = MY_CXT.randcxt;
/*
* Fisher-Yates shuffle with first 'k' selections returned.
*
* There is only one algorithm here, no shortcuts other than
* detecting an empty list.
*
* With a list input, the input is on the stack ST(1),ST(2),...
* We move the last item to ST(0) then shuffle 'k' iterations.
*
* With an array reference input, we cannot modify the input at all.
* We create an index array and shuffle using that. Remembering to
* act like the last item is at the front so we match the list results.
* We optimize by pushing each selection onto the return stack as
* we find it rather than pushing them all at the end with another loop.
*/
if (items > 2 || !SvROK(ST(1)) || SvTYPE(SvRV(ST(1))) != SVt_PVAV) {
/* Standard form, where we are given an array of items */
nitems = items-1;
if (k > nitems)
k = nitems;
if (GIMME_V != G_ARRAY)
XSRETURN_IV(k);
ST(0) = ST(items-1); /* Move last value to the first stack entry. */
for (i = 0; i < k; i++) {
uint32_t j = urandomm32(randcxt, nitems-i);
{ SV* t = ST(i); ST(i) = ST(i+j); ST(i+j) = t; }
}
} else { /* We are given a single array reference. Select from it. */
DECL_ARREF(avp);
USE_ARREF(avp, ST(1), SUBNAME, AR_READ);
nitems = len_avp;
if (nitems == 0)
RETURN_NOTHING();
if (k > nitems)
k = nitems;
if (GIMME_V != G_ARRAY)
XSRETURN_IV(k);
if (nitems < 65536) {
uint16_t *I;
New(0, I, nitems, uint16_t);
I[0] = nitems-1; for (i = 1; i < nitems; i++) I[i] = i-1;
EXTEND(SP, (EXTEND_TYPE)k);
for (i = 0; i < k; i++) {
uint32_t j = urandomm32(randcxt, nitems-i);
uint16_t t = I[i+j]; I[i+j] = I[i];
PUSHs(FETCH_ARREF(avp,t));
}
Safefree(I);
} else {
size_t *I;
New(0, I, nitems, size_t);
I[0] = nitems-1; for (i = 1; i < nitems; i++) I[i] = i-1;
EXTEND(SP, (EXTEND_TYPE)k);
for (i = 0; i < k; i++) {
size_t j = urandomm64(randcxt, nitems-i);
size_t t = I[i+j]; I[i+j] = I[i];
PUSHs(FETCH_ARREF(avp,t));
}
Safefree(I);
}
}
XSRETURN(k);
void is_happy(SV* svn, UV base = 10, UV k = 2)
PREINIT:
UV n, sum;
int h, status;
PPCODE:
if (base < 2 || base > 36) croak("%s: invalid base: %"UVuf, SUBNAME, base);
if (k > 10) croak("%s: invalid exponent %"UVuf, SUBNAME, k);
status = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
if (status == 0 && base == 10) { /* String op to reduce into range. */
STRLEN i, len;
const char* s = SvPV(svn, len);
if (len <= UV_MAX/ipow(9,k)) {
for (sum = 0, i = 0; i < len; i++)
sum += ipow(s[i]-'0',k);
h = happy_height(sum, base, k);
if (h >= 0)
RETURN_NPARITY( (h>0) ? h+1 : 0);
}
}
if (status != 0) {
h = happy_height(n, base, k);
if (h >= 0)
RETURN_NPARITY(h);
}
DISPATCHPP_RETURN();
void
sumdigits(SV* svn, SV* svbase = 0)
PREINIT:
UV sum, base = 10;
STRLEN i, len;
const char* s;
PPCODE:
SvGETMAGIC(svn);
if (!SvOK(svn)) croak("Parameter must be defined");
if (items > 1 && !_validate_and_set(&base, aTHX_ svbase, IFLAG_NONNEG))
DISPATCHPP_RETURN();
if (base < 2 || base > 36) croak("%s: invalid base: %"UVuf, SUBNAME, base);
sum = 0;
/* faster for integer input in base 10 */
if (base == 10 && SVNUMTEST(svn) && (SvIsUV(svn) || SvIVX(svn) >= 0)) {
UV n, t = my_svuv(svn);
while ((n=t)) {
t = n / base;
sum += n - base*t;
}
XSRETURN_UV(sum);
}
s = SvPV(svn, len);
/* If no base given and input is 0x... or 0b..., select base. */
if (items < 2 && len > 2 && s[0] == '0' && (s[1] == 'x' || s[1] == 'b')){
base = (s[1] == 'x') ? 16 : 2;
s += 2;
len -= 2;
}
for (i = 0; i < len; i++) {
UV d = 0;
const char c = s[i];
if (c >= '0' && c <= '9') { d = c - '0'; }
else if (c >= 'a' && c <= 'z') { d = c - 'a' + 10; }
else if (c >= 'A' && c <= 'Z') { d = c - 'A' + 10; }
if (d < base)
sum += d;
}
XSRETURN_UV(sum);
void reverse_digits(SV* svn, SV* svbase = 0)
PREINIT:
int nstatus, bstatus = 1;
UV n, base = 10, r = 0;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
if (items > 1)
bstatus = _validate_and_set(&base, aTHX_ svbase, IFLAG_NONNEG);
if (base < 2) croak("%s: invalid base: %"UVuf, SUBNAME, base);
if (nstatus == 1 && bstatus == 1) {
while (n > 0) {
UV q = n / base;
UV d = n - q * base;
if (r > (UV_MAX - d) / base)
break; /* overflow, dispatch */
r = r * base + d;
n = q;
}
if (n == 0)
XSRETURN_UV(r);
}
if (bstatus == 1) {
char *rstr;
STRLEN len, rlen;
const char *s = SvPV_nomg(svn, len);
rstr = strint_reverse_digits(s, len, base, &rlen);
if (rstr)
RETURN_SIGN_STRINT_STR(1, rstr, rlen);
}
DISPATCHPP_RETURN();
void todigits(SV* svn, SV* svbase = 0, SV* svtlen = 0)
ALIAS:
todigitstring = 1
PREINIT:
int nstatus, bstatus, lstatus, tlen, i;
UV n, base;
PPCODE:
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
bstatus = lstatus = 1;
base = 10;
tlen = -1;
if (items > 1) {
bstatus = _validate_and_set(&base, aTHX_ svbase, IFLAG_NONNEG);
if (base < 2 || (ix == 1 && base > 36))
croak("%s: invalid base: %"UVuf, SUBNAME, base);
}
if (items > 2) {
UV uvtlen;
lstatus = _validate_and_set(&uvtlen, aTHX_ svtlen, IFLAG_NONNEG);
if (lstatus != 0 && uvtlen > (UV)PERL_INT_MAX) lstatus = 0;
if (lstatus != 0) tlen = uvtlen;
}
if (bstatus != 1 || lstatus != 1)
DISPATCHPP_RETURN_GMPIF(0);
if (nstatus != 0 && tlen <= 128) {
if (ix == 0) { /* todigits with native input */
UV digits[128];
int len = to_digit_array(digits, n, base, tlen);
if (len >= 0) {
if (GIMME_V != G_ARRAY) XSRETURN_UV(len);
EXTEND(SP, (EXTEND_TYPE)len);
for (i = 0; i < len; i++)
PUSH_NPARITY( digits[len-i-1] );
XSRETURN(len);
}
} else { /* todigitstring with native input */
char s[128+1];
IV len = to_digit_string(s, n, base, tlen);
if (len >= 0) {
XPUSHs(sv_2mortal(newSVpv(s, len)));
XSRETURN(1);
}
}
}
/* todigits or todigitstring base 10 (large size) */
if (base == 10 && tlen < 0) {
STRLEN len;
char *str = SvPV(svn, len);
if (len > 0 && (*str == '-' || *str == '+')) {
str++;
len--;
}
while (len > 1 && *str == '0') {
str++;
len--;
}
if (len == 1 && str[0] == '0') {
if (ix == 1) {
XPUSHs(sv_2mortal(newSVpvn("", 0)));
XSRETURN(1);
}
RETURN_NOTHING();
}
if (ix == 1) {
XPUSHs(sv_2mortal(newSVpv(str, len)));
XSRETURN(1);
}
if (GIMME_V != G_ARRAY) XSRETURN_UV(len);
EXTEND(SP, (EXTEND_TYPE)len);
for (i = 0; i < (int)len; i++)
PUSH_NPARITY(str[i]-'0');
XSRETURN(len);
}
#if BITS_PER_WORD == 64
DISPATCHPP_RETURN_GMPIF(base <= UINT32_MAX);
#else
DISPATCHPP_RETURN();
#endif
void fromdigits(SV* svn, SV* svbase = 0)
PREINIT:
int bstatus;
UV n, base;
STRLEN slen, rlen;
char *s, *str;
int status;
PPCODE:
SvGETMAGIC(svn);
if (!SvOK(svn)) croak("Parameter must be defined");
bstatus = 1;
base = 10;
if (items > 1)
bstatus = _validate_and_set(&base, aTHX_ svbase, IFLAG_NONNEG);
if (base < 2) croak("%s: invalid base: %"UVuf, SUBNAME, base);
if (bstatus == 1) {
if (SvROK(svn) && SvTYPE(SvRV(svn)) == SVt_PVAV) {
size_t len;
UV* r = 0;
if (arrayref_to_digit_array(aTHX_ &len, &r, svn, base)) {
if (from_digit_to_UV(&n, r, len, base)) {
Safefree(r);
XSRETURN_UV(n);
} else if ((str = strint_fromdigits(r, len, base, &rlen)) != 0) {
Safefree(r);
RETURN_SIGN_STRINT_STR(1, str, rlen);
}
Safefree(r);
}
} else if (!SvROK(svn) || _sv_is_math_object(aTHX_ svn)) {
s = SvPV_nomg(svn, slen);
status = strint_fromdigitstring(&n, &str, &rlen, s, slen, base);
if (status == 1)
XSRETURN_UV(n);
if (status == 2)
RETURN_SIGN_STRINT_STR(1, str, rlen);
croak("fromdigits: internal error");
} else {
croak("fromdigits: first argument must be a string or array reference");
}
}
DISPATCHPP_RETURN();
void is_harshad(SV* svn, SV* svbase = 0)
PREINIT:
int nstatus, bstatus;
UV n, base;
PPCODE:
if (items == 1) { bstatus = 1; base = 10; }
else { bstatus = _validate_and_set(&base, aTHX_ svbase, IFLAG_NONNEG); }
if (bstatus == 1) {
if (base < 2) croak("%s: invalid base: %"UVuf, SUBNAME, base);
if (base <= UINT32_MAX) {
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_ANY);
if (nstatus == -1 || (nstatus == 1 && n == 0))
RETURN_NPARITY(0);
if (nstatus == 1) {
UV N, t, sum;
uint32_t b = (uint32_t) base;
for (sum = 0, N = n; N > 0; N = t) {
t = N / b;
sum += N - b*t;
}
RETURN_NPARITY(n % sum == 0);
}
}
}
/* We can read the string and sum the digits, but no way to mod here. */
DISPATCHPP_RETURN();
void is_palindrome(SV* svn, SV* svbase = 0)
PREINIT:
int bstatus;
UV n, base;
PPCODE:
if (items == 1) { bstatus = 1; base = 10; }
else { bstatus = _validate_and_set(&base, aTHX_ svbase, IFLAG_NONNEG); }
if (bstatus == 1) {
if (base < 2) croak("%s: invalid base: %"UVuf, SUBNAME, base);
if (base <= UINT32_MAX &&
_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
uint32_t b = (uint32_t) base;
UV forward = n, reverse = 0;
while (n > 0) {
uint32_t digit = n % b;
reverse = reverse * b + digit;
n /= b;
}
RETURN_NPARITY(forward == reverse);
}
}
DISPATCHPP_RETURN();
void digital_root(SV* svn, IN SV* svbase = 0)
ALIAS:
mult_digital_root = 1
PREINIT:
UV n, dr, base;
int bstatus;
PPCODE:
if (items==1){bstatus = 1; base = 10;}
else {bstatus = _validate_and_set(&base,aTHX_ svbase,IFLAG_NONNEG);}
if (bstatus == 1 && _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
if (base < 2) croak("%s: invalid base: %"UVuf, SUBNAME, base);
if (n == 0) {
dr = 0;
} else if (ix == 0) {
dr = 1 + (n-1) % (base-1);
} else {
UV digits[128];
int i, len;
dr = n;
while (dr >= base) {
len = to_digit_array(digits, dr, base, -1);
for (dr = 1, i = 0; i < len; i++)
dr *= digits[i];
}
}
RETURN_NPARITY(dr);
}
DISPATCHPP_RETURN();
void tozeckendorf(SV* svn)
PREINIT:
UV n;
PPCODE:
if (_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG)) {
char *str = to_zeckendorf(n);
XPUSHs(sv_2mortal(newSVpv(str, 0)));
Safefree(str);
XSRETURN(1);
}
DISPATCHPP_RETURN();
void fromzeckendorf(IN SV* svstr)
PREINIT:
int status;
STRLEN len;
const char* str;
PPCODE:
SvGETMAGIC(svstr);
if (!SvOK(svstr)) croak("Parameter must be defined");
str = SvPV_nomg(svstr, len);
status = validate_zeckendorf(str, len);
if (status == 0)
croak("fromzeckendorf: expected binary string");
if (status == -1)
croak("fromzeckendorf: expected binary string in canonical Zeckendorf form");
if (status == 1)
XSRETURN_UV(from_zeckendorf(str, len));
DISPATCHPP_RETURN();
void
lastfor()
PREINIT:
dMY_CXT;
PPCODE:
/* printf("last for with count = %u\n", MY_CXT.forcount); */
if (MY_CXT.forcount == 0) croak("lastfor called outside a loop");
MY_CXT.forexit = 1;
/* In some ideal world this would also act like a last */
return;
#define START_FORCOUNT \
do { \
New(0, forguard, 1, forcount_guard_t); \
forguard->previous_forcount = MY_CXT.forcount; \
forguard->expected_forcount = ++MY_CXT.forcount; \
forguard->previous_forexit = MY_CXT.forexit; \
forguard->active = 1; \
forexit = &MY_CXT.forexit; \
*forexit = 0; \
SAVEDESTRUCTOR_X(forcount_guard_cleanup, forguard); \
} while(0)
#define CHECK_FORCOUNT \
if (*forexit) break;
#define END_FORCOUNT \
do { \
current_forcount = MY_CXT.forcount; \
/* Put back outer loop's exit request, if any. */ \
*forexit = forguard->previous_forexit; \
MY_CXT.forcount = forguard->previous_forcount; \
forguard->active = 0; \
/* Ensure loops are nested and not woven. */ \
if (current_forcount != forguard->expected_forcount) croak("for loop mismatch"); \
} while (0)
#define DECL_FORCOUNT \
forcount_guard_t *forguard; \
uint16_t current_forcount; \
char *forexit
void
forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
PROTOTYPE: &$;$
PREINIT:
SV* svarg;
CV *subcv;
unsigned char* segment;
UV beg, end, seg_base, seg_low, seg_high;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&beg, aTHX_ svbeg, IFLAG_NONNEG) ||
(svend && !_validate_and_set(&end, aTHX_ svend, IFLAG_NONNEG)))
DISPATCHPP_RETURN_VOID();
if (!svend) { end = beg; beg = 2; }
START_FORCOUNT;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(beg);
GvSV(PL_defgv) = svarg;
/* Handle early part */
#if USE_MULTICALL
if (!CvISXSUB(subcv) && beg <= end) {
SC_dMULTICALL;
I32 gimme = G_VOID;
SC_PUSH_MULTICALL(subcv);
if (beg < 6) {
beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
for ( ; beg < 6 && beg <= end; beg += 1+(beg>2) ) {
CHECK_FORCOUNT;
sv_setuv(svarg, beg);
SC_MULTICALL;
}
}
if (beg <= end) {
if (
#if BITS_PER_WORD == 64
(beg >= UVCONST( 100000000000000) && end-beg < 100000) ||
(beg >= UVCONST( 10000000000000) && end-beg < 40000) ||
(beg >= UVCONST( 1000000000000) && end-beg < 17000) ||
#endif
((end-beg) < 500) ) { /* MULTICALL next prime */
for (beg = next_prime(beg-1); beg <= end && beg != 0; beg = next_prime(beg)) {
CHECK_FORCOUNT;
sv_setuv(svarg, beg);
SC_MULTICALL;
}
} else { /* MULTICALL segment sieve */
void* ctx = start_segment_primes(beg, end, &segment);
while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
int crossuv = (seg_high > IV_MAX) && !SvIsUV(svarg);
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_base, seg_low, seg_high )
CHECK_FORCOUNT;
/* sv_setuv(svarg, p); */
if (SvTYPE(svarg) != SVt_IV) { sv_setuv(svarg, p); }
else if (crossuv && p > IV_MAX) { sv_setuv(svarg, p); crossuv=0; }
else { SvUV_set(svarg, p); }
SC_MULTICALL;
END_DO_FOR_EACH_SIEVE_PRIME
CHECK_FORCOUNT;
}
end_segment_primes(ctx);
}
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
{
if (beg < 6) {
beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
for ( ; beg < 6 && beg <= end; beg += 1+(beg>2) ) {
sv_setuv(svarg, beg);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
}
}
if (beg <= end) { /* NO-MULTICALL segment sieve */
void* ctx = start_segment_primes(beg, end, &segment);
while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_base, seg_low, seg_high )
CHECK_FORCOUNT;
sv_setuv(svarg, p);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
END_DO_FOR_EACH_SIEVE_PRIME
CHECK_FORCOUNT;
}
end_segment_primes(ctx);
}
}
SvREFCNT_dec(svarg);
END_FORCOUNT;
#define FORCOMPTEST(ix,n) \
( (ix==1) || (ix==0 && n&1) )
void
foroddcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
ALIAS:
forcomposites = 1
PROTOTYPE: &$;$
PREINIT:
UV beg, end;
SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
CV *subcv;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&beg, aTHX_ svbeg, IFLAG_NONNEG) ||
(svend && !_validate_and_set(&end, aTHX_ svend, IFLAG_NONNEG)))
DISPATCHPP_RETURN_VOID();
if (!svend) { end = beg; beg = ix ? 4 : 9; }
START_FORCOUNT;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
GvSV(PL_defgv) = svarg;
#if USE_MULTICALL
if (!CvISXSUB(subcv) && end >= beg) {
unsigned char* segment;
UV seg_base, seg_low, seg_high, c, cbeg, cend, cinc, prevprime, nextprime;
void* ctx;
SC_dMULTICALL;
I32 gimme = G_VOID;
SC_PUSH_MULTICALL(subcv);
if (beg >= MPU_MAX_PRIME ||
#if BITS_PER_WORD == 64
(beg >= UVCONST( 100000000000000) && end-beg < 120000) ||
(beg >= UVCONST( 10000000000000) && end-beg < 50000) ||
(beg >= UVCONST( 1000000000000) && end-beg < 20000) ||
#endif
end-beg < 1000 ) {
beg = (beg <= 4) ? 3 : beg-1;
nextprime = next_prime(beg);
while (beg++ < end) {
if (beg == nextprime)
nextprime = next_prime(beg);
else if (FORCOMPTEST(ix,beg)) {
sv_setuv(svarg, beg);
SC_MULTICALL;
}
CHECK_FORCOUNT;
}
} else {
if (!ix) {
if (beg < 8) beg = 8;
} else if (beg <= 4) { /* sieve starts at 7, so handle this here */
sv_setuv(svarg, 4);
SC_MULTICALL;
beg = 6;
}
/* Find the two primes that bound their interval. */
/* beg must be < max_prime, and end >= max_prime is special. */
prevprime = prev_prime(beg);
nextprime = (end >= MPU_MAX_PRIME) ? MPU_MAX_PRIME : next_prime(end);
ctx = start_segment_primes(beg, nextprime, &segment);
while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
int crossuv = (seg_high > IV_MAX) && !SvIsUV(svarg);
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_base, seg_low, seg_high )
cbeg = prevprime+1;
if (cbeg < beg)
cbeg = beg - (ix == 0 && (beg % 2));
prevprime = p;
cend = prevprime-1; if (cend > end) cend = end;
/* If ix=0, skip evens by starting 1 farther and skipping by 2 */
cinc = 1 + (ix==0);
for (c = cbeg + (ix==0); c <= cend; c += cinc) {
CHECK_FORCOUNT;
if (SvTYPE(svarg) != SVt_IV) { sv_setuv(svarg,c); }
else if (crossuv && c > IV_MAX) { sv_setuv(svarg,c); crossuv=0;}
else { SvUV_set(svarg,c); }
SC_MULTICALL;
}
END_DO_FOR_EACH_SIEVE_PRIME
}
end_segment_primes(ctx);
if (end > nextprime) /* Complete the case where end > max_prime */
while (nextprime++ < end)
if (FORCOMPTEST(ix,nextprime)) {
CHECK_FORCOUNT;
sv_setuv(svarg, nextprime);
SC_MULTICALL;
}
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
if (beg <= end) {
beg = (beg <= 4) ? 3 : beg-1;
while (beg++ < end) {
if (FORCOMPTEST(ix,beg) && !is_prob_prime(beg)) {
sv_setuv(svarg, beg);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
}
}
}
SvREFCNT_dec(svarg);
END_FORCOUNT;
void
forsemiprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
PROTOTYPE: &$;$
PREINIT:
UV beg, end;
SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
CV *subcv;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&beg, aTHX_ svbeg, IFLAG_NONNEG) ||
(svend && !_validate_and_set(&end, aTHX_ svend, IFLAG_NONNEG)))
DISPATCHPP_RETURN_VOID();
if (!svend) { end = beg; beg = 4; }
if (beg < 4) beg = 4;
if (end > MPU_MAX_SEMI_PRIME) end = MPU_MAX_SEMI_PRIME;
START_FORCOUNT;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
GvSV(PL_defgv) = svarg;
#if USE_MULTICALL
if (!CvISXSUB(subcv) && end >= beg) {
UV c, seg_beg, seg_end, *S, count;
SC_dMULTICALL;
I32 gimme = G_VOID;
SC_PUSH_MULTICALL(subcv);
if (beg >= MPU_MAX_SEMI_PRIME ||
#if BITS_PER_WORD == 64
(beg >= UVCONST(10000000000000000000) && end-beg < 1400000) ||
(beg >= UVCONST( 1000000000000000000) && end-beg < 950000) ||
(beg >= UVCONST( 100000000000000000) && end-beg < 440000) ||
(beg >= UVCONST( 10000000000000000) && end-beg < 240000) ||
(beg >= UVCONST( 1000000000000000) && end-beg < 65000) ||
(beg >= UVCONST( 100000000000000) && end-beg < 29000) ||
(beg >= UVCONST( 10000000000000) && end-beg < 11000) ||
(beg >= UVCONST( 1000000000000) && end-beg < 5000) ||
#endif
end-beg < 200 ) {
for (c = beg; c <= end && c >= beg; c++) {
if (is_semiprime(c)) {
sv_setuv(svarg, c);
SC_MULTICALL;
}
CHECK_FORCOUNT;
}
} else {
while (beg < end) {
seg_beg = beg;
seg_end = end;
if ((seg_end - seg_beg) > 50000000) seg_end = seg_beg + 50000000 - 1;
count = range_semiprime_sieve(&S, seg_beg, seg_end);
for (c = 0; c < count; c++) {
sv_setuv(svarg, S[c]);
SC_MULTICALL;
CHECK_FORCOUNT;
}
Safefree(S);
beg = seg_end+1;
CHECK_FORCOUNT;
}
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
if (beg <= end) {
beg = (beg <= 4) ? 3 : beg-1;
while (beg++ < end) {
if (is_semiprime(beg)) {
sv_setuv(svarg, beg);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
}
}
}
SvREFCNT_dec(svarg);
END_FORCOUNT;
void
foralmostprimes (SV* block, IN UV k, IN SV* svbeg, IN SV* svend = 0)
PROTOTYPE: &$$;$
PREINIT:
UV c, beg, end, shiftres;
SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
CV *subcv;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&beg, aTHX_ svbeg, IFLAG_NONNEG) ||
(svend && !_validate_and_set(&end, aTHX_ svend, IFLAG_NONNEG)))
DISPATCHPP_RETURN_VOID();
if (!svend) { end = beg; beg = 1; }
/* If k is over 63 but the beg/end points are UVs, then we're empty. */
if (k == 0 || k >= BITS_PER_WORD) XSRETURN(0);
if (beg < (UVCONST(1) << k)) beg = UVCONST(1) << k;
if (end > max_nth_almost_prime(k)) end = max_nth_almost_prime(k);
if (beg > end) XSRETURN(0);
/* We might be able to reduce the k value. */
shiftres = 0;
if (k > MPU_MAX_POW3)
shiftres = k - MPU_MAX_POW3;
while ((k-shiftres) > 1 && (end >> shiftres) < ipow(3, k - shiftres))
shiftres++;
beg = (beg >> shiftres) + (((beg >> shiftres) << shiftres) < beg);
end = end >> shiftres;
k -= shiftres;
/* k <= 40 (64-bit) or 20 (32-bit). */
START_FORCOUNT;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
GvSV(PL_defgv) = svarg;
#if USE_MULTICALL
if (!CvISXSUB(subcv) && end >= beg) {
UV seg_beg, seg_end, *S, count, k3 = ipow(3,k);
SC_dMULTICALL;
I32 gimme = G_VOID;
SC_PUSH_MULTICALL(subcv);
while (beg <= end) {
/* TODO: Tuning this better would be nice */
UV ssize = 65536 * 256;
seg_beg = beg;
seg_end = end;
if (k > 12) ssize *= 16;
if (k > 18 || seg_beg > 9*k3) ssize *= 4;
if (k > 24 || seg_beg > 81*k3) ssize *= 3;
if ((seg_end - seg_beg) > ssize) seg_end = seg_beg + ssize - 1;
count = generate_almost_primes(&S, k, seg_beg, seg_end);
for (c = 0; c < count; c++) {
sv_setuv(svarg, S[c] << shiftres);
SC_MULTICALL;
CHECK_FORCOUNT;
}
Safefree(S);
if (seg_end == UV_MAX) break;
beg = seg_end+1;
CHECK_FORCOUNT;
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
if (beg <= end) {
for (c = beg; c <= end && c >= beg; c++) {
if (is_almost_prime(k,c)) {
sv_setuv(svarg, c << shiftres);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
}
}
}
SvREFCNT_dec(svarg);
END_FORCOUNT;
void
fordivisors (SV* block, IN SV* svn)
PROTOTYPE: &$
PREINIT:
UV i, n, ndivisors;
UV *divs;
SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
CV *subcv;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG))
DISPATCHPP_RETURN_VOID();
divs = divisor_list(n, &ndivisors, UV_MAX);
START_FORCOUNT;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
GvSV(PL_defgv) = svarg;
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
I32 gimme = G_VOID;
SC_PUSH_MULTICALL(subcv);
for (i = 0; i < ndivisors; i++) {
sv_setuv(svarg, divs[i]);
SC_MULTICALL;
CHECK_FORCOUNT;
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
{
for (i = 0; i < ndivisors; i++) {
sv_setuv(svarg, divs[i]);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
}
}
SvREFCNT_dec(svarg);
Safefree(divs);
END_FORCOUNT;
void
forpart (SV* block, IN SV* svn, IN SV* svh = 0)
ALIAS:
forcomp = 1
PROTOTYPE: &$;$
PREINIT:
UV i, n, amin, amax, nmin, nmax, tmp;
int primeq;
bool doloop0, doloopn;
CV *subcv;
SV** svals;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG))
DISPATCHPP_RETURN_VOID();
if (n > (UV_MAX-2)) croak("%s: argument overflow", SUBNAME);
amin = 1; amax = n; nmin = 1; nmax = n; primeq = -1;
if (svh != 0) {
HV* rhash;
SV** svp;
int argstatus = 1;
if (!SvROK(svh) || SvTYPE(SvRV(svh)) != SVt_PVHV)
croak("%s: expected hash reference", SUBNAME);
rhash = (HV*) SvRV(svh);
if ((svp = hv_fetchs(rhash, "n", 0)) != NULL) {
argstatus &= (_validate_and_set(&nmin, aTHX_ *svp, IFLAG_NONNEG) == 1);
nmax = nmin;
}
if ((svp = hv_fetchs(rhash, "amin", 0)) != NULL)
argstatus &= (_validate_and_set(&amin, aTHX_ *svp, IFLAG_NONNEG) == 1);
if ((svp = hv_fetchs(rhash, "amax", 0)) != NULL)
argstatus &= (_validate_and_set(&amax, aTHX_ *svp, IFLAG_NONNEG) == 1);
if ((svp = hv_fetchs(rhash, "nmin", 0)) != NULL)
argstatus &= (_validate_and_set(&nmin, aTHX_ *svp, IFLAG_NONNEG) == 1);
if ((svp = hv_fetchs(rhash, "nmax", 0)) != NULL)
argstatus &= (_validate_and_set(&nmax, aTHX_ *svp, IFLAG_NONNEG) == 1);
if ((svp = hv_fetchs(rhash, "prime",0)) != NULL) {
argstatus &= (_validate_and_set(&tmp, aTHX_ *svp, IFLAG_NONNEG) == 1);
primeq = (tmp != 0); /* primeq: -1=any, 0=!prime, 1=prime */
}
if (argstatus == 0)
DISPATCHPP_RETURN_VOID();
if (amin < 1) amin = 1;
if (amax > n) amax = n;
if (nmin < 1) nmin = 1;
if (nmax > n) nmax = n;
}
if (primeq == 1) {
UV prev = prev_prime(amax+1);
UV next = amin <= 2 ? 2 : next_prime(amin-1);
if (amin < next) amin = next;
if (amax > prev) amax = prev;
}
doloop0 = (n == 0 && nmin <= 1);
doloopn = (n >= nmin && nmin <= nmax && amin <= amax && nmax>0 && amax>0);
if (!doloop0 && !doloopn) XSRETURN(0);
START_FORCOUNT;
if (doloop0) { /* Empty, no setup needed */
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
}
if (doloopn && !*forexit) {
/* RuleAsc algorithm from Kelleher and O'Sullivan 2009/2014) */
a[k++] = x;
x = r;
y -= x;
}
a[k] = x + y;
/* ------ length restrictions ------ */
while (k+1 > nmax) { /* Skip range if over max size */
a[k-1] += a[k];
k--;
}
/* Look into: quick skip over nmin range */
if (k+1 < nmin) { /* Skip if not over min size */
if (a[0] >= n-nmin+1 && a[k] > 1) break; /* early exit check */
continue;
}
/* ------ value restrictions ------ */
if (amin > 1 || amax < n) {
/* Lexical order allows us to start at amin, and exit early */
if (a[0] > amax) break;
if (ix == 0) { /* value restrictions for partitions */
if (a[k] > amax) continue;
} else { /* restrictions for compositions */
/* TODO: maybe skip forward? */
for (i = 0; i <= k; i++)
if (a[i] < amin || a[i] > amax)
break;
if (i <= k) continue;
}
}
if (primeq != -1) {
for (i = 0; i <= k; i++) if (is_prime(a[i]) != primeq) break;
if (i <= k) continue;
}
PUSHMARK(SP); EXTEND(SP, (EXTEND_TYPE)k+1);
for (i = 0; i <= k; i++) { PUSHs(svals[a[i]]); }
PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
}
Safefree(a);
for (i = 0; i <= n; i++)
SvREFCNT_dec(svals[i]);
Safefree(svals);
}
END_FORCOUNT;
void
forcomb (SV* block, IN SV* svn, IN SV* svk = 0)
PROTOTYPE: &$;$
PREINIT:
UV i, n, k, begk, endk;
int nstatus, kstatus;
CV *subcv;
SV** svals;
UV* cm;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
nstatus = _validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG);
kstatus = 1;
if (nstatus != 0) {
if (items == 2) {
begk = 0; endk = n;
} else if (_validate_and_set(&k, aTHX_ svk, IFLAG_NONNEG)) {
begk = endk = k;
if (begk > n)
XSRETURN(0);
} else {
kstatus = 0;
}
}
if (nstatus == 0 || kstatus == 0)
DISPATCHPP_RETURN_VOID();
New(0, svals, n, SV*);
for (i = 0; i < n; i++) {
svals[i] = newSVuv(i);
SvREADONLY_on(svals[i]);
}
New(0, cm, endk+1, UV);
START_FORCOUNT;
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
I32 gimme = G_VOID;
AV *av = save_ary(PL_defgv);
AvREAL_off(av);
SC_PUSH_MULTICALL(subcv);
for (k = begk; k <= endk; k++) {
_comb_init(cm, k, 0);
while (1) {
IV j;
av_extend(av, k-1);
av_fill(av, k-1);
for (j = k-1; j >= 0; j--)
AvARRAY(av)[j] = svals[ cm[k-j-1]-1 ];
SC_MULTICALL;
CHECK_FORCOUNT;
if (_comb_iterate(cm, k, n, 0)) break;
}
CHECK_FORCOUNT;
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
} else
#endif
{
for (k = begk; k <= endk; k++) {
_comb_init(cm, k, 0);
while (1) {
PUSHMARK(SP); EXTEND(SP, (EXTEND_TYPE)k);
for (i = 0; i < k; i++) { PUSHs(svals[ cm[k-i-1]-1 ]); }
PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
if (_comb_iterate(cm, k, n, 0)) break;
}
CHECK_FORCOUNT;
}
}
Safefree(cm);
for (i = 0; i < n; i++)
SvREFCNT_dec(svals[i]);
Safefree(svals);
END_FORCOUNT;
void forperm (SV* block, IN SV* svn)
ALIAS:
forderange = 1
PROTOTYPE: &$
PREINIT:
UV i, n;
CV *subcv;
SV** svals;
UV* cm;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&n, aTHX_ svn, IFLAG_NONNEG))
DISPATCHPP_RETURN_VOID();
if (n <= 1) {
if (!(ix == 1 && n == 1)) {
START_FORCOUNT;
PUSHMARK(SP); EXTEND(SP, 1);
if (n == 1) PUSH_NPARITY(0);
PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
END_FORCOUNT;
}
XSRETURN(0);
}
New(0, svals, n, SV*);
for (i = 0; i < n; i++) {
svals[i] = newSVuv(i);
SvREADONLY_on(svals[i]);
}
New(0, cm, n+1, UV);
START_FORCOUNT;
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
I32 gimme = G_VOID;
AV *av = save_ary(PL_defgv);
AvREAL_off(av);
SC_PUSH_MULTICALL(subcv);
_comb_init(cm, n, ix == 1);
while (1) {
IV j;
av_extend(av, n-1);
av_fill(av, n-1);
for (j = n-1; j >= 0; j--)
AvARRAY(av)[j] = svals[ cm[n-j-1]-1 ];
SC_MULTICALL;
CHECK_FORCOUNT;
if (_comb_iterate(cm, n, n, 1+ix)) break;
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
} else
#endif
{
_comb_init(cm, n, ix == 1);
while (1) {
PUSHMARK(SP); EXTEND(SP, (EXTEND_TYPE)n);
for (i = 0; i < n; i++) { PUSHs(svals[ cm[n-i-1]-1 ]); }
PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
if (_comb_iterate(cm, n, n, 1+ix)) break;
}
}
Safefree(cm);
for (i = 0; i < n; i++)
SvREFCNT_dec(svals[i]);
Safefree(svals);
END_FORCOUNT;
void forsetproduct (SV* block, ...)
PROTOTYPE: &@
PREINIT:
SSize_t narrays, i, j, *arlen, *arcnt;
SV ***arsvs;
CV *subcv;
forsetproduct_guard_t *fsguard;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
narrays = items-1;
if (narrays < 1) XSRETURN(0);
for (i = 1; i <= narrays; i++) {
SvGETMAGIC(ST(i));
CHECK_ARRAYREF(ST(i));
if (av_count((AV *)SvRV(ST(i))) == 0)
XSRETURN(0);
}
Newz(0, fsguard, 1, forsetproduct_guard_t);
fsguard->narrays = narrays;
fsguard->active = 1;
SAVEDESTRUCTOR_X(forsetproduct_guard_cleanup, fsguard);
Newz(0, arcnt, narrays, SSize_t); fsguard->arcnt = arcnt;
Newz(0, arlen, narrays, SSize_t); fsguard->arlen = arlen;
Newz(0, arsvs, narrays, SV**); fsguard->arsvs = arsvs;
/* Make local SV copies. Allows magic/tied inputs and prevents source aliasing. */
for (i = 0; i < narrays; i++) {
DECL_ARREF(inav);
USE_ARREF(inav, ST(i+1), SUBNAME, AR_READ);
arlen[i] = len_inav;
Newz(0, arsvs[i], len_inav, SV*);
for (j = 0; j < (SSize_t)len_inav; j++) {
SV* v = FETCH_ARREF(inav,j);
arsvs[i][j] = v ? newSVsv(v) : newSV(0);
REFRESH_ARREF(inav);
}
}
START_FORCOUNT;
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
SV **arr;
I32 gimme = G_VOID;
AV *av = save_ary(PL_defgv);
SSize_t lastarg = narrays-1;
AvREAL_off(av);
SC_PUSH_MULTICALL(subcv);
av_extend(av, lastarg);
av_fill(av, lastarg);
arr = AvARRAY(av);
while (1) {
for (i = lastarg; i >= 0; i--)
arr[i] = arsvs[i][arcnt[i]];
SC_MULTICALL;
CHECK_FORCOUNT;
for (i = lastarg; i >= 0; i--) {
if (++arcnt[i] >= arlen[i]) arcnt[i] = 0;
else break;
}
if (i < 0)
break;
if (av_count(av) != (Size_t)narrays || AvARRAY(av) != arr) {
av_fill(av, lastarg);
arr = AvARRAY(av);
}
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
do {
PUSHMARK(SP); EXTEND(SP, (EXTEND_TYPE)narrays);
for (i = 0; i < narrays; i++) { PUSHs(arsvs[i][arcnt[i]]); }
PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
CHECK_FORCOUNT;
for (i = narrays-1; i >= 0; i--) {
if (++arcnt[i] >= arlen[i]) arcnt[i] = 0;
else break;
}
} while (i >= 0);
forsetproduct_guard_release(aTHX_ fsguard);
END_FORCOUNT;
void
forfactored (SV* block, IN SV* svbeg, IN SV* svend = 0)
ALIAS:
forsquarefree = 1
PROTOTYPE: &$;$
PREINIT:
UV beg, end, n, *factors;
int i, nfactors, maxfactors;
factor_range_context_t fctx;
SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
CV *subcv;
SV* svals[64];
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&beg, aTHX_ svbeg, IFLAG_NONNEG) ||
(svend && !_validate_and_set(&end, aTHX_ svend, IFLAG_NONNEG)))
DISPATCHPP_RETURN_VOID();
if (!svend) { end = beg; beg = 1; }
if (beg < 1) beg = 1;
if (beg > end) XSRETURN(0);
START_FORCOUNT;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
GvSV(PL_defgv) = svarg;
if (beg <= 1) {
sv_setuv(svarg, 1);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
beg = 2;
}
if (*forexit || beg > end) {
SvREFCNT_dec(svarg);
END_FORCOUNT;
XSRETURN(0);
}
for (maxfactors = 0, n = end >> 1; n; n >>= 1)
maxfactors++;
for (i = 0; i < maxfactors; i++) {
svals[i] = newSVuv(UV_MAX);
SvREADONLY_on(svals[i]);
}
fctx = factor_range_init(beg, end, ix);
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
I32 gimme = G_VOID;
AV *av = save_ary(PL_defgv);
AvREAL_off(av);
SC_PUSH_MULTICALL(subcv);
for (n = 0; n < end-beg+1; n++) {
CHECK_FORCOUNT;
nfactors = factor_range_next(&fctx);
if (nfactors > 0) {
sv_setuv(svarg, fctx.n);
factors = fctx.factors;
av_extend(av, nfactors-1);
av_fill(av, nfactors-1);
for (i = nfactors-1; i >= 0; i--) {
SV* sv = svals[i];
SvREADONLY_off(sv);
sv_setuv(sv, factors[i]);
SvREADONLY_on(sv);
AvARRAY(av)[i] = sv;
}
SC_MULTICALL;
}
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
for (n = 0; n < end-beg+1; n++) {
CHECK_FORCOUNT;
nfactors = factor_range_next(&fctx);
if (nfactors > 0) {
PUSHMARK(SP); EXTEND(SP, (EXTEND_TYPE)nfactors);
sv_setuv(svarg, fctx.n);
factors = fctx.factors;
for (i = 0; i < nfactors; i++) {
SV* sv = svals[i];
SvREADONLY_off(sv);
sv_setuv(sv, factors[i]);
SvREADONLY_on(sv);
PUSHs(sv);
}
PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
}
}
factor_range_destroy(&fctx);
SvREFCNT_dec(svarg);
for (i = 0; i < maxfactors; i++)
SvREFCNT_dec(svals[i]);
END_FORCOUNT;
void forsquarefreeint(SV* block, IN SV* svbeg, IN SV* svend = 0)
PROTOTYPE: &$;$
PREINIT:
UV beg, end, i;
unsigned char* isf;
SV* svarg; /* We use svarg to prevent clobbering $_ outside the block */
CV *subcv;
DECL_FORCOUNT;
dMY_CXT;
PPCODE:
SETSUBREF(subcv, block);
if (!_validate_and_set(&beg, aTHX_ svbeg, IFLAG_NONNEG) ||
(svend && !_validate_and_set(&end, aTHX_ svend, IFLAG_NONNEG)))
DISPATCHPP_RETURN_VOID();
if (!svend) { end = beg; beg = 1; }
if (beg < 1) beg = 1;
if (beg > end) XSRETURN(0);
START_FORCOUNT;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
GvSV(PL_defgv) = svarg;
if (beg <= 1) {
sv_setuv(svarg, 1);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
beg = 2;
}
if (*forexit) {
SvREFCNT_dec(svarg);
END_FORCOUNT;
XSRETURN(0);
}
while (beg <= end) {
UV seglo = beg, seghi = end;
if (seghi-seglo > (65536*256))
seghi = seglo + 65536*256 - 1;
isf = range_issquarefree(seglo, seghi);
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
I32 gimme = G_VOID;
SC_PUSH_MULTICALL(subcv);
for (i = 0; i < seghi-seglo+1; i++) {
CHECK_FORCOUNT;
if (isf[i]) {
sv_setuv(svarg, seglo+i);
SC_MULTICALL;
}
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
for (i = 0; i < seghi-seglo+1; i++) {
CHECK_FORCOUNT;
if (isf[i]) {
sv_setuv(svarg, seglo+i);
PUSHMARK(SP); PUTBACK; call_sv((SV*)subcv, G_VOID|G_DISCARD); SPAGAIN;
}
}
Safefree(isf);
if (seghi == UV_MAX) break;
beg = seghi+1;
CHECK_FORCOUNT;
}
SvREFCNT_dec(svarg);
END_FORCOUNT;
SV *atmp, *btmp;
SETSUBREF(subcv, block);
if (items <= 2)
RETURN_NOTHING();
New(0, retsvarr, items-2, SV*);
atmp = sv_newmortal();
btmp = sv_newmortal();
agv = gv_fetchpv("a", GV_ADD, SVt_PV);
bgv = gv_fetchpv("b", GV_ADD, SVt_PV);
SAVESPTR(GvSV(agv));
SAVESPTR(GvSV(bgv));
GvSV(agv) = atmp;
GvSV(bgv) = btmp;
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
SV *cbret;
I32 gimme = G_SCALAR;
SC_PUSH_MULTICALL(subcv);
for (i = 1; i < items-1; i++) {
SvSetMagicSV(atmp, args[i]);
SvSetMagicSV(btmp, args[i+1]);
SC_MULTICALL_SCALAR(cbret);
retsvarr[i-1] = newSVsv(cbret);
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
{
for (i = 1; i < items-1; i++) {
SvSetMagicSV(atmp, args[i]);
SvSetMagicSV(btmp, args[i+1]);
retsvarr[i-1] = newSVsv(xs_call_cv_noinput_1_sv(aTHX_ subcv));
}
}
if (GIMME_V != G_ARRAY) {
for (i = 0; i < items-2; i++)
{ SvREFCNT_dec_NN(retsvarr[i]); retsvarr[i]=0; }
Safefree(retsvarr);
XSRETURN_UV(items-2);
}
for (i = 0; i < items-2; i++)
{ ST(i) = sv_2mortal(retsvarr[i]); retsvarr[i]=0; }
Safefree(retsvarr);
XSRETURN(items-2);
}
void
vecwindow(SV* block, SV* svstep, SV* svsize, ...)
PROTOTYPE: &$$@
PREINIT:
SSize_t i, j, nlist, step, size, numret;
UV ustep, usize;
CV *subcv;
SV **list;
AV *list_av, *result_av;
PPCODE:
{ /* Generalized sliding window: block, step, size, list.
Each window of 'size' elements is passed as @_ to block;
all return values are collected and returned (like map). */
SETSUBREF(subcv, block);
if (!_validate_and_set(&ustep, aTHX_ svstep, IFLAG_POS) ||
!_validate_and_set(&usize, aTHX_ svsize, IFLAG_POS) ||
ustep > (UV)MAX_SSIZET || usize > (UV)MAX_SSIZET) {
DISPATCHPP_RETURN();
}
step = (SSize_t)ustep;
size = (SSize_t)usize;
/* ST(0)=block, ST(1)=step, ST(2)=size, ST(3..)=list */
nlist = items - 3;
if (nlist < size)
RETURN_NOTHING();
/* Save the input SV* -- callbacks could realloc or reuse the arg stack. */
list_av = (AV*)sv_2mortal((SV*)newAV());
av_extend(list_av, nlist-1);
for (i = 0; i < nlist; i++)
av_push(list_av, SvREFCNT_inc(PL_stack_base[ax+3+i]));
list = AvARRAY(list_av);
/* result_av is where we're going to push the callback results. */
result_av = (AV*)sv_2mortal((SV*)newAV());
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
I32 gimme = G_ARRAY;
AV *av = save_ary(PL_defgv);
SC_PUSH_MULTICALL(subcv);
for (i = 0; i + size <= nlist; i += step) {
av_clear(av);
for (j = 0; j < size; j++)
av_push(av, newSVsv(list[i + j]));
SC_MULTICALL_ARRAY(result_av);
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
{
for (i = 0; i + size <= nlist; i += step) {
I32 k, nret;
PUSHMARK(SP); EXTEND(SP, (EXTEND_TYPE)size);
for (j = 0; j < size; j++) PUSHs(sv_2mortal(newSVsv(list[i + j])));
PUTBACK;
nret = call_sv((SV*)subcv, G_ARRAY);
SPAGAIN;
for (k = 0; k < nret; k++)
av_push(result_av, newSVsv(SP[k + 1 - nret]));
SP -= nret; PUTBACK;
}
}
numret = (SSize_t)av_count(result_av);
if (GIMME_V == G_ARRAY) {
SV **res = AvARRAY(result_av);
AvREAL_off(result_av); /* transfer ownership to stack mortals */
EXTEND(SP, (EXTEND_TYPE)numret);
for (i = 0; i < numret; i++)
PUSHs(sv_2mortal(res[i]));
XSRETURN(numret);
}
XSRETURN_UV(numret);
}
void
vecpairwise(SV* block, SV* sva, SV* svb)
PROTOTYPE: &$$
PREINIT:
SSize_t i, npairs, numret;
CV *subcv;
AV *result_av;
GV *agv, *bgv;
SV *atmp, *btmp;
DECL_ARREF(arr1);
DECL_ARREF(arr2);
PPCODE:
{ /* Similar to pairwise from List::MoreUtils, but block values do not alias inputs. */
SETSUBREF(subcv, block);
if (!SvROK(sva) || SvTYPE(SvRV(sva)) != SVt_PVAV ||
!SvROK(svb) || SvTYPE(SvRV(svb)) != SVt_PVAV)
croak("vecpairwise: expected two array references");
USE_ARREF(arr1, sva, SUBNAME, AR_READ);
USE_ARREF(arr2, svb, SUBNAME, AR_READ);
npairs = (len_arr1 < len_arr2) ? (SSize_t)len_arr1 : (SSize_t)len_arr2;
if (npairs <= 0)
RETURN_NOTHING();
result_av = (AV*)sv_2mortal((SV*)newAV());
atmp = sv_newmortal();
btmp = sv_newmortal();
agv = gv_fetchpv("a", GV_ADD, SVt_PV);
bgv = gv_fetchpv("b", GV_ADD, SVt_PV);
SAVESPTR(GvSV(agv));
SAVESPTR(GvSV(bgv));
GvSV(agv) = atmp;
GvSV(bgv) = btmp;
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
I32 gimme = G_ARRAY;
SC_PUSH_MULTICALL(subcv);
for (i = 0; i < npairs; i++) {
SV *a = FETCH_ARREF(arr1, i);
SV *b = FETCH_ARREF(arr2, i);
SvSetMagicSV(atmp, a);
SvSetMagicSV(btmp, b);
SC_MULTICALL_ARRAY(result_av);
REFRESH_ARREF(arr1);
REFRESH_ARREF(arr2);
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
{
for (i = 0; i < npairs; i++) {
I32 k, nret;
SV *a = FETCH_ARREF(arr1, i);
SV *b = FETCH_ARREF(arr2, i);
SvSetMagicSV(atmp, a);
SvSetMagicSV(btmp, b);
PUSHMARK(SP); PUTBACK;
nret = call_sv((SV*)subcv, G_ARRAY);
SPAGAIN;
for (k = 0; k < nret; k++)
av_push(result_av, newSVsv(SP[k + 1 - nret]));
SP -= nret; PUTBACK;
REFRESH_ARREF(arr1);
REFRESH_ARREF(arr2);
}
}
numret = (SSize_t)av_count(result_av);
if (GIMME_V == G_ARRAY) {
SV **res = AvARRAY(result_av);
AvREAL_off(result_av); /* transfer ownership to stack mortals */
EXTEND(SP, (EXTEND_TYPE)numret);
for (i = 0; i < numret; i++)
PUSHs(sv_2mortal(res[i]));
XSRETURN(numret);
}
XSRETURN_UV(numret);
}
void
vecnone(SV* block, ...)
ALIAS:
vecall = 1
vecany = 2
vecnotall = 3
vecfirst = 4
vecfirstidx = 6
PROTOTYPE: &@
PPCODE:
{ /* This is very similar to List::Util. Try to maintain compat. */
int ret_true = !(ix & 2); /* return true at end of loop for none/all; false for any/notall */
int invert = (ix & 1); /* invert block test for all/notall */
SSize_t index;
SV **args = &PL_stack_base[ax];
CV *subcv;
SETSUBREF(subcv, block);
SAVESPTR(GvSV(PL_defgv));
#if USE_MULTICALL
if (!CvISXSUB(subcv)) {
SC_dMULTICALL;
SV *cbret;
I32 gimme = G_SCALAR;
SC_PUSH_MULTICALL(subcv);
for (index = 1; index < items; index++) {
GvSV(PL_defgv) = args[index];
SC_MULTICALL_SCALAR(cbret);
if (SvTRUEx(cbret) ^ invert)
break;
}
FIX_MULTICALL_REFCOUNT;
SC_POP_MULTICALL;
}
else
#endif
{
for (index = 1; index < items; index++) {
GvSV(PL_defgv) = args[index];
if (SvTRUEx(xs_call_cv_noinput_1_sv(aTHX_ subcv)) ^ invert)
break;
}
}
if (ix == 4) {
if (index == items)
XSRETURN_UNDEF;
ST(0) = ST(index);
XSRETURN(1);
}
if (ix == 6) {
if (index == items)
XSRETURN_IV(-1);
XSRETURN_UV(index-1);
}
if (index != items) /* We exited the loop early */
ret_true = !ret_true;
RETURN_NPARITY(ret_true ? 1 : 0);
}
void vecuniq(...)
PROTOTYPE: @
PREINIT:
int typemask, retvals;
SSize_t j;
PPCODE:
retvals = (GIMME_V != G_SCALAR && GIMME_V != G_VOID);
/* Check all inputs to see if we can quickly process as an iset. */
for (j = 0, typemask = 0; j < items; j++) {
SV *sv = ST(j);
SvGETMAGIC(sv);
if (!SvOK(sv) || !SVNUMTEST(sv))
break;
if ((UV)SvIVX(sv) > (UV)IV_MAX)
typemask |= (SvIsUV(sv) ? ISET_TYPE_UV : ISET_TYPE_IV);
if (typemask == ISET_TYPE_INVALID)
break;
}
/* 2. Use an iset if possible. */
if (j >= items) {
iset_t s;
size_t sz, nret = 0;
iset_create(&s, (size_t)items);
for (j = 0, nret = 0; j < items; j++) {
SV *sv = ST(j);
UV n = (UV)SvIVX(sv);
int nsign = (SvIsUV(sv) || SvIVX(sv) >= 0) ? 1 : -1;
if (iset_add(&s, n, nsign) && retvals) {
PUSHs(ST(j));
nret++;
}
}
if (iset_sign(&s) == 0) {
iset_destroy(&s);
croak("vecuniq: iset failure");
}
sz = iset_size(&s);
iset_destroy(&s);
if (!retvals)
XSRETURN_UV(sz);
if (nret != sz)
croak("vecuniq: iset %lu items, pushed %lu items", (unsigned long)sz, (unsigned long)nret);
XSRETURN(nret);
}
/* 3. The generic Perl hash method */
{
/* Based on List::MoreUtils::XS by Parseval and Rehsack */
UV count = 0;
HV *hv = newHV();
SV **args = &PL_stack_base[ax];
SV *tmp = sv_newmortal();
sv_2mortal(newRV_noinc((SV*)hv));
for (j = 0; j < items; j++) {
SvGETMAGIC(args[j]);
if (!SvOK(args[j]))
croak("vecuniq: all values must be defined");
if (retvals) SvSetSV_nosteal(tmp, args[j]);
else sv_setsv_nomg(tmp, args[j]);
if (!hv_exists_ent(hv, tmp, 0)) {
if (retvals) args[count] = args[j];
count++;
hv_store_ent(hv, tmp, &PL_sv_yes, 0);
}
}
if (retvals)
XSRETURN(count);
XSRETURN_UV(count);
}
void vecfreq(...)
PROTOTYPE: @
PREINIT:
int itype;
size_t len, i, retlen;
UV *L, count;
PPCODE:
if (items == 0)
RETURN_NOTHING();
/* Try to read native integers. Bail to PP if something else. */
len = (size_t) items;
New(0, L, len, UV);
itype = IARR_TYPE_ANY;
for (i = 0; i < len && itype != IARR_TYPE_BAD && SVNUMTEST(ST(i)); i++) {
IV n = SvIVX(ST(i));
if (n < 0) {
if (SvIsUV(ST(i))) itype |= IARR_TYPE_POS;
else itype |= IARR_TYPE_NEG;
}
L[i] = n;
}
if (i < len || itype == IARR_TYPE_BAD) {
Safefree(L);
DISPATCHPP_RETURN();
}
if (itype == IARR_TYPE_NEG)
sort_iv_array((IV*)L, len);
else
sort_uv_array(L, len);
/* 2. Walk the sorted integers */
if (GIMME_V == G_SCALAR) {
count = 0;
for (i = 1; i < len; i++)
if (L[i] != L[i-1])
count++;
ST(0) = sv_2mortal(newSVuv(count+1));
retlen = 1;
} else {
int sign = itype == IARR_TYPE_NEG ? -1 : 1;
EXTEND(SP, (EXTEND_TYPE)len*2);
retlen = 0;
count = 1;
for (i = 1; i < len; i++) {
if (L[i] == L[i-1]) { count++; continue; }
PUSHs(sv_2mortal(NEWSVINT(sign,L[i-1]))); /* key */
PUSHs(sv_2mortal(newSVuv(count))); /* val */
retlen += 2;
count = 1;
}
PUSHs(sv_2mortal(NEWSVINT(sign,L[i-1]))); /* key */
PUSHs(sv_2mortal(newSVuv(count))); /* val */
retlen += 2;
}
Safefree(L);
XSRETURN(retlen);
void vecsingleton(...)
PROTOTYPE: @
PREINIT:
int itype;
size_t len, i, retlen, count;
UV *L;
iset_t seen, dups;
PPCODE:
if (items == 0)
RETURN_NOTHING();
/* Try to read native integers. Bail to PP if something else. */
len = (size_t) items;
New(0, L, len, UV);
iset_create(&seen, len);
iset_create(&dups, len>>1);
itype = IARR_TYPE_ANY;
for (i = 0; i < len && itype != IARR_TYPE_BAD && SVNUMTEST(ST(i)); i++) {
IV n = SvIVX(ST(i));
int sign = 1;
if (n < 0) {
if (SvIsUV(ST(i))) itype |= IARR_TYPE_POS;
else { itype |= IARR_TYPE_NEG; sign = -1; }
}
L[i] = n;
if (!iset_add(&seen, n, sign))
iset_add(&dups, n, sign);
}
if (iset_is_invalid(&seen)) itype = IARR_TYPE_BAD; /* Poison the type */
iset_destroy(&seen);
if (i < len || itype == IARR_TYPE_BAD) {
iset_destroy(&dups);
Safefree(L);
DISPATCHPP_RETURN();
}
if (GIMME_V != G_ARRAY) {
for (i = 0, count = 0; i < len; i++)
if (!iset_contains(&dups, L[i]))
count++;
ST(0) = sv_2mortal(newSVuv(count));
retlen = 1;
} else {
for (i = 0, retlen = 0; i < len; i++)
if (!iset_contains(&dups, L[i]))
ST(retlen++) = ST(i);
}
iset_destroy(&dups);
Safefree(L);
XSRETURN(retlen);
( run in 2.071 seconds using v1.01-cache-2.11-cpan-4e7a2411597 )