[libmath-prime-util-perl] 48/59: Factoring speedups
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:45:02 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.10
in repository libmath-prime-util-perl.
commit 2d1025e0c9c882bbd311de205fc0017f4e85d3d4
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Jul 13 02:44:03 2012 -0600
Factoring speedups
---
XS.xs | 69 +++++++++++++++++++------------------
factor.c | 118 ++++++++++++++++++++++++++-------------------------------------
2 files changed, 84 insertions(+), 103 deletions(-)
diff --git a/XS.xs b/XS.xs
index fb0085b..7413c56 100644
--- a/XS.xs
+++ b/XS.xs
@@ -199,6 +199,7 @@ _XS_factor(IN UV n)
if (n < 4) {
XPUSHs(sv_2mortal(newSVuv( n ))); /* If n is 0-3, we're done. */
} else {
+ int const verbose = 0;
UV tlim = 53; /* Below this we've checked with trial division */
UV tofac_stack[MPU_MAX_FACTORS+1];
UV factored_stack[MPU_MAX_FACTORS+1];
@@ -227,46 +228,48 @@ _XS_factor(IN UV n)
do { /* loop over each remaining factor */
while ( (n >= (tlim*tlim)) && (!is_definitely_prime(n)) ) {
int split_success = 0;
- /* For larger n, do a different sequence */
- if (n > UVCONST(10000000) ) { /* tune this */
- /* SQUFOF (succeeds 98-99.9%) */
- if (!split_success) {
- split_success = squfof_factor(n, tofac_stack+ntofac, 256*1024)-1;
-//printf("squfof %d\n", split_success);
- }
- /* A few rounds of Pollard's Rho usually gets the factors */
- if (!split_success) {
- split_success = prho_factor(n, tofac_stack+ntofac, 800)-1;
-//printf("prho %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;
-//printf("holf %d\n", split_success);
- }
- /* Less than 0.00003% of numbers make it past here. */
- if (!split_success) {
- split_success = prho(n, tofac_stack+ntofac, 256*1024)-1;
-//printf("prho %d\n", split_success);
- }
- } else {
- /* A few rounds of Pollard rho is good for finding small factors */
- if (!split_success) {
- split_success = prho_factor(n, tofac_stack+ntofac, 800)-1;
-//if (split_success) printf("small prho 1: %llu %llu\n", tofac_stack[ntofac], tofac_stack[ntofac+1]); else printf("small prho 0\n");
- }
+
+ if (!split_success) {
+ split_success = pbrent_factor(n, tofac_stack+ntofac, 1500)-1;
+ if (verbose) { if (split_success) printf("pbrent 1: %lu %lu\n", tofac_stack[ntofac], tofac_stack[ntofac+1]); else printf("pbrent 0\n"); }
+ }
+ if (!split_success) {
+ split_success = squfof_factor(n, tofac_stack+ntofac, 256*1024)-1;
+ if (verbose) printf("squfof %d\n", split_success);
+ }
+
+ /* Time to do expensive primality check and exit now if it is */
+ if (!split_success && _XS_is_prime(n)) {
+ if (verbose) printf("oops, %lu is prime\n", n);
+ factored_stack[nfactored++] = n;
+ n = 1;
+ break;
+ }
+
+ /* A few rounds of Pollard's Rho usually gets the factors */
+ if (!split_success) {
+ split_success = prho_factor(n, tofac_stack+ntofac, 800)-1;
+ if (verbose) printf("prho %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);
}
+
+ /* A miniscule fraction of numbers 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);
+ }
+
if (split_success) {
MPUassert( split_success == 1, "split factor returned more than 2 factors");
ntofac++; /* Leave one on the to-be-factored stack */
n = tofac_stack[ntofac]; /* Set n to the other one */
- } else if (_XS_is_prime(n)) {
- /* No wonder we couldn't factor it. It's prime. */
- factored_stack[nfactored++] = n;
- n = 1;
} else {
/* trial divisions */
-printf("doing trial on %llu\n", n);
+ if (verbose) printf("doing trial on %llu\n", n);
UV f = tlim;
UV m = tlim % 30;
UV limit = (UV) (sqrt(n)+0.1);
diff --git a/factor.c b/factor.c
index e8621cd..f6be9d6 100644
--- a/factor.c
+++ b/factor.c
@@ -397,7 +397,7 @@ int fermat_factor(UV n, UV *factors, UV rounds)
}
/* Hart's One Line Factorization.
- * Translation from my GMP, missing perfect square calc and premult.
+ * Missing premult (hard to do in native precision without overflow)
*/
int holf_factor(UV n, UV *factors, UV rounds)
{
@@ -407,7 +407,9 @@ int holf_factor(UV n, UV *factors, UV rounds)
for (i = 1; i <= rounds; i++) {
s = sqrt( (double)n * (double)i );
- if ( (s*s) != (n*i) ) s++;
+ /* Assume s^2 isn't a perfect square. We're rapidly losing precision
+ * so we won't be able to accurately detect it anyway. */
+ s++; /* s = ceil(sqrt(n*i)) */
m = sqrmod(s, n);
if (is_perfect_square(m, &f)) {
f = gcd_ui( (s>f) ? s-f : f-s, n);
@@ -425,94 +427,71 @@ int holf_factor(UV n, UV *factors, UV rounds)
}
-/* Pollard / Brent
- *
- * Probabilistic. If you give this a prime number, it will loop
- * until it runs out of rounds.
- */
+/* Pollard / Brent. Brent's modifications to Pollard's Rho. Maybe faster. */
int pbrent_factor(UV n, UV *factors, UV rounds)
{
- UV a, f, i;
+ UV f, i, r;
UV Xi = 2;
UV Xm = 2;
+ UV a = 1;
+ const UV inner = 128;
MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in pbrent_factor");
- switch (n%4) {
- case 0: a = 1; break;
- case 1: a = 3; break;
- case 2: a = 5; break;
- case 3: a = 7; break;
- default: a = 11; break;
- }
-
- for (i = 1; i <= rounds; i++) {
- Xi = sqraddmod(Xi, a, n);
- f = gcd_ui( (Xi>Xm) ? Xi-Xm : Xm-Xi, n);
- if ( (f != 1) && (f != n) ) {
- factors[0] = f;
- factors[1] = n/f;
- MPUassert( factors[0] * factors[1] == n , "incorrect factoring");
- return 2;
+ r = 1;
+ while (rounds > 0) {
+ UV rleft = (r > rounds) ? rounds : r;
+ UV saveXi;
+ /* Do rleft rounds, inner at a time */
+ while (rleft > 0) {
+ UV dorounds = (rleft > inner) ? inner : rleft;
+ UV m = 1;
+ saveXi = Xi;
+ for (i = 0; i < dorounds; i++) {
+ Xi = sqraddmod(Xi, a, n);
+ f = (Xi>Xm) ? Xi-Xm : Xm-Xi;
+ m = mulmod(m, f, n);
+ }
+ rleft -= dorounds;
+ rounds -= dorounds;
+ f = gcd_ui(m, n);
+ if (f != 1)
+ break;
}
- if ( (i & (i-1)) == 0) /* i is a power of 2 */
+ /* If f == 1, then we didn't find a factor. Move on. */
+ if (f == 1) {
+ r *= 2;
Xm = Xi;
- }
- factors[0] = n;
- return 1;
-}
-
-/* Pollard's Rho
- *
- * Probabilistic. If you give this a prime number, it will loop
- * until it runs out of rounds.
- */
-#if 0
-int prho_factor(UV n, UV *factors, UV rounds)
-{
- int in_loop = 0;
- UV a, f, i;
- UV U = 7;
- UV V = 7;
-
- MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in prho_factor");
-
- switch (n%4) {
- case 0: a = 5; break;
- case 1: a = 7; break;
- case 2: a = 11; break;
- case 3: a = 1; break;
- default: a = 3; break;
- }
-
- for (i = 1; i <= rounds; i++) {
- U = sqraddmod(U, a, n);
- V = sqraddmod(V, a, n);
- V = sqraddmod(V, a, n);
- f = gcd_ui( (U > V) ? U-V : V-U, n);
- if (f == n) {
- if (in_loop++) /* Mark that we've been here */
- break; /* Exit now if we're cycling */
- } else if (f != 1) {
- factors[0] = f;
- factors[1] = n/f;
- MPUassert( factors[0] * factors[1] == n , "incorrect factoring");
- return 2;
+ continue;
}
+ if (f == n) { /* backup, with safety */
+ Xi = saveXi;
+ do {
+ Xi = sqraddmod(Xi, a, n);
+ f = gcd_ui( (Xi>Xm) ? Xi-Xm : Xm-Xi, n);
+ } while (f == 1 && r-- != 0);
+ if ( (f == 1) || (f == n) ) break;
+ }
+ factors[0] = f;
+ factors[1] = n/f;
+ MPUassert( factors[0] * factors[1] == n , "incorrect factoring");
+ return 2;
}
factors[0] = n;
return 1;
}
-#else
+
+/* Pollard's Rho. */
int prho_factor(UV n, UV *factors, UV rounds)
{
UV a, f, i, m, oldU, oldV;
- const UV inner = 8;
+ const UV inner = 128;
UV U = 7;
UV V = 7;
MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in prho_factor");
+ /* We could just as well say a = 1 */
switch (n%8) {
case 1: a = 1; break;
case 3: a = 2; break;
@@ -535,7 +514,7 @@ int prho_factor(UV n, UV *factors, UV rounds)
f = gcd_ui(m, n);
if (f == 1)
continue;
- if (f == n) { /* backup */
+ if (f == n) { /* back up to find a factor*/
U = oldU; V = oldV;
i = inner;
do {
@@ -555,7 +534,6 @@ int prho_factor(UV n, UV *factors, UV rounds)
factors[0] = n;
return 1;
}
-#endif
/* Pollard's P-1
*
--
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