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 <d...@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 _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits