[libmath-prime-util-perl] 71/181: Move divisor_sum to XS; more XS interface cleanup
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:08 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit 053dfe07cf14e67b004d4ec1e379dbb73f10b8ce
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Dec 30 15:09:25 2013 -0800
Move divisor_sum to XS; more XS interface cleanup
---
XS.xs | 153 +++++++++++++++++++++-------------------------
factor.c | 27 ++++----
factor.h | 2 +-
lib/Math/Prime/Util.pm | 81 ++----------------------
lib/Math/Prime/Util/PP.pm | 69 +++++++++++++++++++++
5 files changed, 159 insertions(+), 173 deletions(-)
diff --git a/XS.xs b/XS.xs
index fb3a9ca..851142f 100644
--- a/XS.xs
+++ b/XS.xs
@@ -38,8 +38,8 @@
#define my_svuv(sv) PSTRTOULL(SvPV_nolen(sv), NULL, 10)
#define my_sviv(sv) PSTRTOLL(SvPV_nolen(sv), NULL, 10)
#else
- #define my_svuv(sv) (!SvROK(sv)) ? SvUV(sv) : PSTRTOULL(SvPV_nolen(sv), NULL, 10)
- #define my_sviv(sv) (!SvROK(sv)) ? SvIV(sv) : PSTRTOLL(SvPV_nolen(sv), NULL, 10)
+ #define my_svuv(sv) ((!SvROK(sv)) ? SvUV(sv) : PSTRTOULL(SvPV_nolen(sv), NULL, 10))
+ #define my_sviv(sv) ((!SvROK(sv)) ? SvIV(sv) : PSTRTOLL(SvPV_nolen(sv), NULL, 10))
#endif
/* multicall compatibility stuff */
@@ -148,42 +148,50 @@ PROTOTYPES: ENABLE
void
-_XS_set_verbose(IN int verbose)
-
-int
-_XS_get_verbose()
-
-void
-_XS_set_callgmp(IN int call_gmp)
-
-int
-_XS_get_callgmp()
-
+prime_memfree()
+ ALIAS:
+ _prime_memfreeall = 1
+ _XS_get_verbose = 2
+ _XS_get_callgmp = 3
+ _get_prime_cache_size = 4
+ _XS_prime_maxbits = 5
+ PPCODE:
+ UV ret = 0;
+ switch (ix) {
+ case 0: prime_memfree(); break;
+ case 1: _prime_memfreeall(); break;
+ case 2: ret = _XS_get_verbose(); break;
+ case 3: ret = _XS_get_callgmp(); break;
+ case 4: ret = get_prime_cache(0,0); break;
+ case 5: ret = BITS_PER_WORD; break;
+ }
+ if (ix > 1) XSRETURN_UV(ret);
void
prime_precalc(IN UV n)
+ ALIAS:
+ _XS_set_verbose = 1
+ _XS_set_callgmp = 2
+ PPCODE:
+ switch (ix) {
+ case 0: prime_precalc(n); break;
+ case 1: _XS_set_verbose(n); break;
+ default: _XS_set_callgmp(n); break;
+ }
-void
-prime_memfree()
void
-_prime_memfreeall()
-
-UV
_XS_prime_count(IN UV low, IN UV high = 0)
- CODE:
+ PPCODE:
if (high == 0) { /* Without a Perl layer in front of this, we'll have */
high = low; /* the pathological case of a-0 turning into 0-a. */
low = 0;
}
if (GIMME_V == G_VOID) {
prime_precalc(high);
- RETVAL = 0;
} else {
- RETVAL = _XS_prime_count(low, high);
+ PUSHs(sv_2mortal(newSVuv( _XS_prime_count(low, high) )));
}
- OUTPUT:
- RETVAL
UV
_XS_nth_prime(IN UV n)
@@ -193,8 +201,8 @@ _XS_nth_prime(IN UV n)
_XS_legendre_pi = 3
_XS_meissel_pi = 4
_XS_lehmer_pi = 5
- _XS_LMO_pi = 6
- _XS_LMOS_pi = 7
+ _XS_LMOS_pi = 6
+ _XS_LMO_pi = 7
PREINIT:
UV ret;
CODE:
@@ -205,36 +213,17 @@ _XS_nth_prime(IN UV n)
case 3: ret = _XS_legendre_pi(n); break;
case 4: ret = _XS_meissel_pi(n); break;
case 5: ret = _XS_lehmer_pi(n); break;
- case 6: ret = _XS_LMO_pi(n); break;
- case 7: ret = _XS_LMOS_pi(n); break;
- default: croak("_XS_nth_prime: Unknown function alias"); break;
+ case 6: ret = _XS_LMOS_pi(n); break;
+ default:ret = _XS_LMO_pi(n); break;
}
RETVAL = ret;
OUTPUT:
RETVAL
UV
-_XS_divisor_sum(IN UV n, IN UV k)
-
-UV
_XS_legendre_phi(IN UV x, IN UV a)
-UV
-_get_prime_cache_size()
- CODE:
- RETVAL = get_prime_cache(0, 0);
- OUTPUT:
- RETVAL
-
-int
-_XS_prime_maxbits()
- CODE:
- RETVAL = BITS_PER_WORD;
- OUTPUT:
- RETVAL
-
-
SV*
sieve_primes(IN UV low, IN UV high)
ALIAS:
@@ -307,7 +296,7 @@ void
_XS_divisors(IN UV n)
PPCODE:
if (GIMME_V == G_SCALAR) {
- PUSHs(sv_2mortal(newSVuv( _XS_divisor_sum(n, 0) )));
+ PUSHs(sv_2mortal(newSVuv( divisor_sum(n, 0) )));
} else {
UV i, ndivisors;
UV* divs = _divisor_list(n, &ndivisors);
@@ -427,8 +416,7 @@ _XS_is_prime(IN UV n)
case 4: ret = _XS_is_lucas_pseudoprime(n, 2); break;
case 5: ret = _XS_is_frobenius_underwood_pseudoprime(n); break;
case 6: ret = _XS_is_aks_prime(n); break;
- case 7: ret = _XS_BPSW(n); break;
- default: croak("_XS_is_prime: Unknown function alias"); break;
+ default:ret = _XS_BPSW(n); break;
}
RETVAL = ret;
OUTPUT:
@@ -463,7 +451,7 @@ is_prime(IN SV* svn)
sub = (ix == 0) ? "Math::Prime::Util::_generic_is_prime"
: "Math::Prime::Util::_generic_is_prob_prime";
_vcallsub(sub);
- XSRETURN(1);
+ return; /* skip implicit PUTBACK */
}
void
@@ -482,7 +470,7 @@ next_prime(IN SV* svn)
}
_vcallsub((ix == 0) ? "Math::Prime::Util::_generic_next_prime" :
"Math::Prime::Util::_generic_prev_prime" );
- XSRETURN(1);
+ return; /* skip implicit PUTBACK */
void
factor(IN SV* svn)
@@ -505,14 +493,27 @@ factor(IN SV* svn)
return; /* skip implicit PUTBACK */
}
+void
+divisor_sum(IN SV* svn, ...)
+ PPCODE:
+ SV* svk = (items > 1) ? ST(1) : 0;
+ int nstatus = _validate_int(aTHX_ svn, 0);
+ int kstatus = (items == 1 || (SvIOK(svk) && SvIV(svk))) ? 1 : 0;
+ if (nstatus == 1 && kstatus == 1) {
+ UV n = my_svuv(svn);
+ UV k = (items > 1) ? my_svuv(svk) : 1;
+ UV sigma = divisor_sum(n, k);
+ if (sigma != 0) XSRETURN_UV(sigma); /* sigma 0 means overflow */
+ }
+ _vcallsubn(aTHX_ G_SCALAR, "Math::Prime::Util::_generic_divisor_sum",items);
+ return; /* skip implicit PUTBACK */
+
void
znorder(IN SV* sva, IN SV* svn)
- PREINIT:
- int astatus, nstatus;
PPCODE:
- astatus = _validate_int(aTHX_ sva, 0);
- nstatus = _validate_int(aTHX_ svn, 0);
+ int astatus = _validate_int(aTHX_ sva, 0);
+ int nstatus = _validate_int(aTHX_ svn, 0);
if (astatus == 1 && nstatus == 1) {
UV a = my_svuv(sva);
UV n = my_svuv(svn);
@@ -525,10 +526,8 @@ znorder(IN SV* sva, IN SV* svn)
void
znprimroot(IN SV* svn)
- PREINIT:
- int status;
PPCODE:
- status = _validate_int(aTHX_ svn, 1);
+ int status = _validate_int(aTHX_ svn, 1);
if (status != 0) {
UV r, n = my_svuv(svn);
if (status == -1)
@@ -538,15 +537,13 @@ znprimroot(IN SV* svn)
XSRETURN_UV(r);
}
_vcallsub("Math::Prime::Util::_generic_znprimroot");
- XSRETURN(1);
+ return; /* skip implicit PUTBACK */
void
kronecker(IN SV* sva, IN SV* svb)
- PREINIT:
- int astatus, bstatus;
PPCODE:
- astatus = _validate_int(aTHX_ sva, 2);
- bstatus = _validate_int(aTHX_ svb, 2);
+ int astatus = _validate_int(aTHX_ sva, 2);
+ int bstatus = _validate_int(aTHX_ svb, 2);
if (astatus == 1 && bstatus == 1) {
UV a = my_svuv(sva);
UV b = my_svuv(svb);
@@ -577,8 +574,7 @@ _XS_ExponentialIntegral(IN SV* x)
case 2: ret = (double) ld_riemann_zeta(SvNV(x)); break;
case 3: ret = _XS_RiemannR(SvNV(x)); break;
case 4: ret = _XS_chebyshev_theta(SvUV(x)); break;
- case 5: ret = _XS_chebyshev_psi(SvUV(x)); break;
- default: ret = 0;
+ default:ret = _XS_chebyshev_psi(SvUV(x)); break;
}
RETVAL = ret;
OUTPUT:
@@ -621,7 +617,7 @@ euler_phi(IN SV* svlo, ...)
Safefree(totients);
}
} else {
- _vcallsubn(aTHX_ G_ARRAY, "Math::Prime::Util::_generic_euler_phi", 2);
+ _vcallsubn(aTHX_ G_ARRAY,"Math::Prime::Util::_generic_euler_phi",items);
return; /* skip implicit PUTBACK */
}
} else {
@@ -659,7 +655,7 @@ moebius(IN SV* svlo, ...)
Safefree(mu);
}
} else {
- _vcallsubn(aTHX_ G_ARRAY, "Math::Prime::Util::_generic_moebius",2);
+ _vcallsubn(aTHX_ G_ARRAY, "Math::Prime::Util::_generic_moebius",items);
return; /* skip implicit PUTBACK */
}
} else {
@@ -687,7 +683,7 @@ carmichael_lambda(IN SV* svn)
}
_vcallsub( (ix == 0) ? "Math::Prime::Util::_generic_carmichael_lambda"
: "Math::Prime::Util::_generic_exp_mangoldt" );
- XSRETURN(1);
+ return; /* skip implicit PUTBACK */
int
_validate_num(SV* svn, ...)
@@ -721,7 +717,6 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
PPCODE:
{
#if !USE_MULTICALL
- dSP;
PUSHMARK(SP);
XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
PUTBACK;
@@ -738,13 +733,12 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
UV seg_base, seg_low, seg_high;
if (!_validate_int(aTHX_ svbeg, 0) || (items >= 3 && !_validate_int(aTHX_ svend,0))) {
- dSP;
PUSHMARK(SP);
XPUSHs(block); XPUSHs(svbeg); if (svend) XPUSHs(svend);
PUTBACK;
(void) call_pv("Math::Prime::Util::_generic_forprimes", G_VOID|G_DISCARD);
SPAGAIN;
- XSRETURN(0);
+ return;
}
if (items < 3) {
beg = 2;
@@ -754,7 +748,7 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
end = my_svuv(svend);
}
if (beg > end)
- XSRETURN(0);
+ return;
cv = sv_2cv(block, &stash, &gv, 0);
if (cv == Nullcv)
@@ -764,7 +758,6 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
GvSV(PL_defgv) = svarg;
/* Handle early part */
while (beg < 6) {
- dSP;
beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
if (beg <= end) {
sv_setuv(svarg, beg);
@@ -807,7 +800,6 @@ forprimes (SV* block, IN SV* svbeg, IN SV* svend = 0)
POP_MULTICALL;
} else {
START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base ) {
- dSP;
sv_setuv(svarg, seg_base + p);
PUSHMARK(SP);
call_sv((SV*)cv, G_VOID|G_DISCARD);
@@ -833,7 +825,7 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
if (cv == Nullcv)
croak("Not a subroutine reference");
- if (items <= 1) XSRETURN(0);
+ if (items <= 1) return;
if (!_validate_int(aTHX_ svbeg, 0) || (items >= 3 && !_validate_int(aTHX_ svend,0))) {
PUSHMARK(SP);
@@ -841,7 +833,7 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
PUTBACK;
(void) call_pv("Math::Prime::Util::_generic_forcomposites", G_VOID|G_DISCARD);
SPAGAIN;
- XSRETURN(0);
+ return;
}
if (items < 3) {
@@ -852,7 +844,7 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
end = my_svuv(svend);
}
if (beg > end)
- XSRETURN(0);
+ return;
SAVESPTR(GvSV(PL_defgv));
svarg = newSVuv(0);
@@ -895,7 +887,6 @@ forcomposites (SV* block, IN SV* svbeg, IN SV* svend = 0)
beg = (beg <= 4) ? 3 : beg-1;
while (beg++ < end) {
if (!_XS_is_prob_prime(beg)) {
- dSP;
sv_setuv(svarg, beg);
PUSHMARK(SP);
call_sv((SV*)cv, G_VOID|G_DISCARD);
@@ -919,16 +910,15 @@ fordivisors (SV* block, IN SV* svn)
if (cv == Nullcv)
croak("Not a subroutine reference");
- if (items <= 1) XSRETURN(0);
+ if (items <= 1) return;
if (!_validate_int(aTHX_ svn, 0)) {
- dSP;
PUSHMARK(SP);
XPUSHs(block); XPUSHs(svn);
PUTBACK;
(void) call_pv("Math::Prime::Util::_generic_fordivisors", G_VOID|G_DISCARD);
SPAGAIN;
- XSRETURN(0);
+ return;
}
n = my_svuv(svn);
@@ -953,7 +943,6 @@ fordivisors (SV* block, IN SV* svn)
#endif
{
for (i = 0; i < ndivisors; i++) {
- dSP;
sv_setuv(svarg, divs[i]);
PUSHMARK(SP);
call_sv((SV*)cv, G_VOID|G_DISCARD);
diff --git a/factor.c b/factor.c
index 9315961..112a151 100644
--- a/factor.c
+++ b/factor.c
@@ -373,18 +373,25 @@ UV* _divisor_list(UV n, UV *num_divisors)
}
-/*
- * The usual method, on OEIS for instance, is:
+/* The usual method, on OEIS for instance, is:
* (p^(k*(e+1))-1) / (p^k-1)
* but that overflows quicky. Instead we rearrange as:
* 1 + p^k + p^k^2 + ... p^k^e
+ * Return 0 if the result overflowed.
*/
-UV _XS_divisor_sum(UV n, UV k)
+static const UV sigma_overflow[5] =
+#if BITS_PER_WORD == 64
+ {3000000000000000000, 3000000000, 2487240, 64260, 7026};
+#else
+ { 845404560, 52560, 1548, 252, 84};
+#endif
+UV divisor_sum(UV n, UV k)
{
UV factors[MPU_MAX_FACTORS+1];
int nfac, i, j;
UV product = 1;
+ if (k > 5 || (k > 0 && n >= sigma_overflow[k-1])) return 0;
if (n == 0) return (k == 0) ? 2 : 1; /* divisors are [0,1] */
if (n == 1) return 1; /* divisors are [1] */
nfac = factor(n, factors);
@@ -784,15 +791,11 @@ int pplus1_factor(UV n, UV *factors, UV B1)
/* My modification of Ben Buhrow's modification of Bob Silverman's SQUFOF code.
*/
-static IV qqueue[100+1];
-static IV qpoint;
-static void enqu(IV q, IV *iter) {
- qqueue[qpoint] = q;
- if (++qpoint >= 100) *iter = -1;
-}
int squfof_factor(UV n, UV *factors, UV rounds)
{
+ IV qqueue[100+1];
+ IV qpoint;
IV rounds2 = (IV) (rounds/16);
UV temp;
IV iq,ll,l2,p,pnext,q,qlast,r,s,t,i;
@@ -827,10 +830,8 @@ int squfof_factor(UV n, UV *factors, UV rounds)
iq = (s + p)/q;
pnext = iq*q - p;
if (q <= ll) {
- if ((q & 1) == 0)
- enqu(q/2,&jter);
- else if (q <= l2)
- enqu(q,&jter);
+ if ((q & 1) == 0) { qqueue[qpoint] = q/2; if (++qpoint>=100) jter = -1; }
+ else if (q <= l2) { qqueue[qpoint] = q; if (++qpoint>=100) jter = -1; }
if (jter < 0) {
factors[0] = n; return 1;
}
diff --git a/factor.h b/factor.h
index 41cf3c5..7ee466a 100644
--- a/factor.h
+++ b/factor.h
@@ -7,6 +7,7 @@
extern int factor(UV n, UV *factors);
extern int factor_exp(UV n, UV *factors, UV* exponents);
+extern UV divisor_sum(UV n, UV k);
extern int trial_factor(UV n, UV *factors, UV maxtrial);
@@ -19,7 +20,6 @@ extern int pplus1_factor(UV n, UV *factors, UV B);
extern int squfof_factor(UV n, UV *factors, UV rounds);
extern int racing_squfof_factor(UV n, UV *factors, UV rounds);
-extern UV _XS_divisor_sum(UV n, UV k);
extern UV* _divisor_list(UV n, UV *num_divisors);
#endif
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 0fc5705..28f5a19 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -107,6 +107,7 @@ BEGIN {
*moebius = \&Math::Prime::Util::_generic_moebius;
*carmichael_lambda = \&Math::Prime::Util::_generic_carmichael_lambda;
*kronecker = \&Math::Prime::Util::_generic_kronecker;
+ *divisor_sum = \&Math::Prime::Util::_generic_divisor_sum;
*znorder = \&Math::Prime::Util::_generic_znorder;
*znprimroot = \&Math::Prime::Util::_generic_znprimroot;
*factor = \&Math::Prime::Util::_generic_factor;
@@ -1350,84 +1351,10 @@ sub jordan_totient {
return $totient;
}
-my @_ds_overflow = # We'll use BigInt math if the input is larger than this.
- (~0 > 4294967295)
- ? (124, 3000000000000000000, 3000000000, 2487240, 64260, 7026)
- : ( 50, 845404560, 52560, 1548, 252, 84);
-sub divisor_sum {
- my($n, $k) = @_;
- return 1 if defined $n && $n == 1;
+sub _generic_divisor_sum {
+ my($n) = @_;
_validate_num($n) || _validate_positive_integer($n);
-
- # Call the XS routine for k=0 and k=1 immediately if possible.
- if ($n <= $_XS_MAXVAL) {
- if (defined $k) {
- return _XS_divisor_sum($n, 0) if $k eq '0';
- return _XS_divisor_sum($n, 1) if $k eq '1' && $n < $_ds_overflow[1];
- } else {
- return _XS_divisor_sum($n, 1) if $n < $_ds_overflow[1];
- }
- }
-
- if (defined $k && ref($k) eq 'CODE') {
- my $sum = $n-$n;
- if (ref($n) eq 'Math::BigInt') {
- # If the original number was a bigint, make sure all divisors are.
- fordivisors { $sum += $k->(Math::BigInt->new("$_")); } $n;
- } else {
- fordivisors { $sum += $k->($_); } $n;
- }
- return $sum;
- }
-
- croak "Second argument must be a code ref or number"
- unless !defined $k || _validate_num($k) || _validate_positive_integer($k);
- $k = 1 if !defined $k;
-
- my $will_overflow = ($k == 0) ? (length($n) >= $_ds_overflow[0])
- : ($k <= 5) ? ($n >= $_ds_overflow[$k])
- : 1;
-
- return _XS_divisor_sum($n, $k) if $n <= $_XS_MAXVAL && !$will_overflow;
-
- # The standard way is:
- # my $pk = $f ** $k; $product *= ($pk ** ($e+1) - 1) / ($pk - 1);
- # But we get less overflow using:
- # my $pk = $f ** $k; $product *= $pk**E for E in 0 .. e
- # Also separate BigInt and do fiddly bits for better performance.
-
- my $product = 1;
- my @pe = ($n <= $_XS_MAXVAL) ? _XS_factor_exp($n) : factor_exp($n);
-
- if (!$will_overflow) {
- foreach my $f (@pe) {
- my ($p, $e) = @$f;
- if ($k == 0) {
- $product *= ($e+1);
- } else {
- my $pk = $p ** $k;
- my $fmult = $pk + 1;
- foreach my $E (2 .. $e) { $fmult += $pk**$E }
- $product *= $fmult;
- }
- }
- } else {
- $product = Math::BigInt->bone;
- my $bik = Math::BigInt->new("$k");
- foreach my $f (@pe) {
- my ($p, $e) = @$f;
- my $pk = Math::BigInt->new("$p")->bpow($bik);
- if ($e == 1) { $pk->binc(); $product->bmul($pk); }
- elsif ($e == 2) { $pk->badd($pk*$pk)->binc(); $product->bmul($pk); }
- else {
- my $fmult = $pk;
- foreach my $E (2 .. $e) { $fmult += $pk->copy->bpow($E) }
- $fmult->binc();
- $product *= $fmult;
- }
- }
- }
- return $product;
+ return Math::Prime::Util::PP::divisor_sum(@_);
}
# Need proto for the block
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 2f4b7c7..4ae17b2 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -546,6 +546,75 @@ sub moebius_range {
return @mu;
}
+my @_ds_overflow = # We'll use BigInt math if the input is larger than this.
+ (~0 > 4294967295)
+ ? (124, 3000000000000000000, 3000000000, 2487240, 64260, 7026)
+ : ( 50, 845404560, 52560, 1548, 252, 84);
+sub divisor_sum {
+ my($n, $k) = @_;
+ return 1 if defined $n && $n == 1;
+
+ if (defined $k && ref($k) eq 'CODE') {
+ my $sum = $n-$n;
+ if (ref($n) eq 'Math::BigInt') {
+ # If the original number was a bigint, make sure all divisors are.
+ foreach my $d (Math::Prime::Util::divisors($n)) {
+ $sum += $k->(Math::BigInt->new("$d"));
+ }
+ } else {
+ foreach my $d (Math::Prime::Util::divisors($n)) {
+ $sum += $k->($d);
+ }
+ }
+ return $sum;
+ }
+
+ croak "Second argument must be a code ref or number"
+ unless !defined $k || _validate_num($k) || _validate_positive_integer($k);
+ $k = 1 if !defined $k;
+
+ my $will_overflow = ($k == 0) ? (length($n) >= $_ds_overflow[0])
+ : ($k <= 5) ? ($n >= $_ds_overflow[$k])
+ : 1;
+
+ # The standard way is:
+ # my $pk = $f ** $k; $product *= ($pk ** ($e+1) - 1) / ($pk - 1);
+ # But we get less overflow using:
+ # my $pk = $f ** $k; $product *= $pk**E for E in 0 .. e
+ # Also separate BigInt and do fiddly bits for better performance.
+
+ my $product = 1;
+ if (!$will_overflow) {
+ foreach my $f (Math::Prime::Util::factor_exp($n)) {
+ my ($p, $e) = @$f;
+ if ($k == 0) {
+ $product *= ($e+1);
+ } else {
+ my $pk = $p ** $k;
+ my $fmult = $pk + 1;
+ foreach my $E (2 .. $e) { $fmult += $pk**$E }
+ $product *= $fmult;
+ }
+ }
+ } else {
+ $product = Math::BigInt->bone;
+ my $bik = Math::BigInt->new("$k");
+ foreach my $f (Math::Prime::Util::factor_exp($n)) {
+ my ($p, $e) = @$f;
+ my $pk = Math::BigInt->new("$p")->bpow($bik);
+ if ($e == 1) { $pk->binc(); $product->bmul($pk); }
+ elsif ($e == 2) { $pk->badd($pk*$pk)->binc(); $product->bmul($pk); }
+ else {
+ my $fmult = $pk;
+ foreach my $E (2 .. $e) { $fmult += $pk->copy->bpow($E) }
+ $fmult->binc();
+ $product *= $fmult;
+ }
+ }
+ }
+ return $product;
+}
+
#############################################################################
# Lehmer prime count
#
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list