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