[libmath-prime-util-perl] 33/50: Changes for p-1 factoring
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:37 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.13
in repository libmath-prime-util-perl.
commit 6f831a105c0dde3fd020fd51a001d2606dffe167
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Nov 8 01:10:05 2012 -0800
Changes for p-1 factoring
---
XS.xs | 40 ++++++++++++++++++++++++---------
factor.c | 77 +++++++++++++++++++++++++++++++++++++++++-----------------------
factor.h | 2 +-
3 files changed, 80 insertions(+), 39 deletions(-)
diff --git a/XS.xs b/XS.xs
index cc2de61..25dadce 100644
--- a/XS.xs
+++ b/XS.xs
@@ -233,14 +233,14 @@ _XS_factor(IN UV n)
UV br_rounds = ((n>>29) < 100000) ? 1500 : 1500;
UV sq_rounds = 80000; /* 20k 91%, 40k 98%, 80k 99.9%, 120k 99.99% */
- /* Small factors will be found quite rapidly with this */
+ /* About 94% of random inputs are factored with this pbrent call */
if (!split_success) {
split_success = pbrent_factor(n, tofac_stack+ntofac, br_rounds)-1;
if (verbose) { if (split_success) printf("pbrent 1: %"UVuf" %"UVuf"\n", tofac_stack[ntofac], tofac_stack[ntofac+1]); else printf("pbrent 0\n"); }
}
if (!split_success && n < (UV_MAX>>3)) {
- /* SQUFOF does very well with what's left after TD and Rho. */
+ /* SQUFOF with these parameters gets 95% of what's left. */
split_success = racing_squfof_factor(n, tofac_stack+ntofac, sq_rounds)-1;
if (verbose) printf("squfof %d\n", split_success);
}
@@ -251,17 +251,19 @@ _XS_factor(IN UV n)
if (verbose) printf("prho %d\n", split_success);
}
+ /* This p-1 gets about 2/3 of what makes it through the above */
+ if (!split_success) {
+ split_success = pminus1_factor(n, tofac_stack+ntofac, 4000, 40000)-1;
+ if (verbose) printf("pminus1 %d\n", split_success);
+ }
+
/* Some rounds of HOLF, good for close to perfect squares */
if (!split_success) {
split_success = holf_factor(n, tofac_stack+ntofac, 2000)-1;
if (verbose) printf("holf %d\n", split_success);
}
- if (!split_success) {
- split_success = pminus1_factor(n, tofac_stack+ntofac, 1000)-1;
- if (verbose) printf("pminus1 %d\n", split_success);
- }
- /* A miniscule fraction of numbers make it here */
+ /* Less than 0.1% of random inputs make it here */
if (!split_success) {
split_success = prho_factor(n, tofac_stack+ntofac, 256*1024)-1;
if (verbose) printf("long prho %d\n", split_success);
@@ -274,7 +276,7 @@ _XS_factor(IN UV n)
croak("bad factor\n");
n = tofac_stack[ntofac]; /* Set n to the other one */
} else {
- /* trial divisions */
+ /* Factor via trial division. Nothing should make it here. */
UV f = tlim;
UV m = tlim % 30;
UV limit = (UV) (sqrt(n)+0.1);
@@ -371,9 +373,27 @@ prho_factor(IN UV n, IN UV maxrounds = 4*1024*1024)
SIMPLE_FACTOR(prho_factor, n, maxrounds);
void
-pminus1_factor(IN UV n, IN UV maxrounds = 1*1024*1024)
+pminus1_factor(IN UV n, IN UV B1 = 1*1024*1024, IN UV B2 = 0)
PPCODE:
- SIMPLE_FACTOR(pminus1_factor, n, maxrounds);
+ if (B2 == 0)
+ B2 = 10*B1;
+ if (n <= 1) {
+ XPUSHs(sv_2mortal(newSVuv( n )));
+ } else {
+ while ( (n% 2) == 0 ) { n /= 2; XPUSHs(sv_2mortal(newSVuv( 2 ))); }
+ while ( (n% 3) == 0 ) { n /= 3; XPUSHs(sv_2mortal(newSVuv( 3 ))); }
+ while ( (n% 5) == 0 ) { n /= 5; XPUSHs(sv_2mortal(newSVuv( 5 ))); }
+ if (n == 1) { /* done */ }
+ else if (_XS_is_prime(n)) { XPUSHs(sv_2mortal(newSVuv( n ))); }
+ else {
+ UV factors[MPU_MAX_FACTORS+1];
+ int i, nfactors;
+ nfactors = pminus1_factor(n, factors, B1, B2);
+ for (i = 0; i < nfactors; i++) {
+ XPUSHs(sv_2mortal(newSVuv( factors[i] )));
+ }
+ }
+ }
int
_XS_miller_rabin(IN UV n, ...)
diff --git a/factor.c b/factor.c
index c20d246..2d3d68e 100644
--- a/factor.c
+++ b/factor.c
@@ -549,21 +549,19 @@ int prho_factor(UV n, UV *factors, UV rounds)
return 1;
}
-/* Pollard's P-1
- *
- * Probabilistic. If you give this a prime number, it will loop
- * until it runs out of rounds.
- */
-int pminus1_factor(UV n, UV *factors, UV B1)
+/* Pollard's P-1 */
+int pminus1_factor(UV n, UV *factors, UV B1, UV B2)
{
- UV f, q, restartq, restarta;
+ UV q, f;
UV a = 2;
+ UV savea = 2;
+ UV saveq = 2;
+ UV j = 1;
UV sqrtB1 = sqrt(B1);
MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in pminus1_factor");
- restartq = 2;
- restarta = a;
+
for (q = 2; q <= sqrtB1; q = _XS_next_prime(q)) {
- UV k = q;
+ UV k = q*q;
UV kmin = B1/q;
while (k <= kmin)
k *= q;
@@ -572,46 +570,70 @@ int pminus1_factor(UV n, UV *factors, UV B1)
if (a == 0) { factors[0] = n; return 1; }
f = gcd_ui(a-1, n);
if (f == 1) {
- restartq = q;
- restarta = a;
- for ( ; q <= B1; q = _XS_next_prime(q)) {
+ savea = a;
+ saveq = q;
+ for (; q <= B1; q = _XS_next_prime(q)) {
a = powmod(a, q, n);
+ if ( (j++ % 32) == 0) {
+ if (a == 0 || gcd_ui(a-1, n) != 1)
+ break;
+ savea = a;
+ saveq = q;
+ }
}
+ if (a == 0) { factors[0] = n; return 1; }
f = gcd_ui(a-1, n);
}
- /* See if we found more than one factor in stage 1, repeat if so */
+ /* If we found more than one factor in stage 1, backup and single step */
if (f == n) {
- a = restarta;
- for (q = restartq; q <= B1; q = _XS_next_prime(q)) {
+ a = savea;
+ for (q = saveq; q <= B1; q = _XS_next_prime(q)) {
UV k = q;
UV kmin = B1/q;
while (k <= kmin)
k *= q;
a = powmod(a, k, n);
f = gcd_ui(a-1, n);
- if (f != 1 && f != n)
+ if (f != 1)
break;
}
+ /* If f == n again, we could do:
+ * for (savea = 3; f == n && savea < 100; savea = _XS_next_prime(savea)) {
+ * a = savea;
+ * for (q = 2; q <= B1; q = _XS_next_prime(q)) {
+ * ...
+ * }
+ * }
+ * but this could be a huge time sink if B1 is large, so just fail.
+ */
}
- if (f == 1 || f == n) { /* stage 2 */
- UV B2 = B1 * 10;
+
+ /* STAGE 2 */
+ if (f == 1 && B2 > B1) {
UV bm = a;
UV b = 1;
- UV j = 1;
+ UV bmdiff;
UV precomp_bm[111] = {0}; /* Enough for B2 = 189M */
/* calculate (a^q)^2, (a^q)^4, etc. */
- precomp_bm[0] = sqrmod(bm, n);
+ bmdiff = sqrmod(bm, n);
+ precomp_bm[0] = bmdiff;
+ for (j = 1; j < 20; j++) {
+ bmdiff = mulmod(bmdiff,bm,n);
+ bmdiff = mulmod(bmdiff,bm,n);
+ precomp_bm[j] = bmdiff;
+ }
a = powmod(a, q, n);
+ j = 1;
while (q <= B2) {
UV lastq = q;
- UV qdiff, bmdiff;
+ UV qdiff;
q = _XS_next_prime(q);
/* compute a^q = a^lastq * a^(q-lastq) */
qdiff = (q - lastq) / 2 - 1;
if (qdiff >= 111) {
- bmdiff = powmod(bm, q-lastq, n); /* Too big of gap */
+ bmdiff = powmod(bm, q-lastq, n); /* Big gap */
} else {
bmdiff = precomp_bm[qdiff];
if (bmdiff == 0) {
@@ -622,13 +644,12 @@ int pminus1_factor(UV n, UV *factors, UV B1)
precomp_bm[qdiff] = bmdiff;
}
}
- a = mulmod( a, bmdiff, n);
- if (a <= 1) break;
- b = mulmod( b, a-1, n);
- if (b == 0) b = 1;
+ a = mulmod(a, bmdiff, n);
+ if (a == 0) break;
+ b = mulmod(b, a-1, n); /* if b == 0, we found multiple factors */
if ( (j++ % 64) == 0 ) {
f = gcd_ui(b, n);
- if (f != 1 && f != n)
+ if (f != 1)
break;
}
}
diff --git a/factor.h b/factor.h
index 1b6f876..d2267fd 100644
--- a/factor.h
+++ b/factor.h
@@ -18,7 +18,7 @@ extern int pbrent_factor(UV n, UV *factors, UV maxrounds);
extern int prho_factor(UV n, UV *factors, UV maxrounds);
-extern int pminus1_factor(UV n, UV *factors, UV maxrounds);
+extern int pminus1_factor(UV n, UV *factors, UV B1, UV B2);
extern int _XS_miller_rabin(UV n, const UV *bases, int nbases);
extern int _XS_is_prob_prime(UV n);
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list