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 8347ffa7a637d4fed148fd17009910246af61da4 Author: Dana Jacobsen <d...@acm.org> Date: Thu Jul 12 18:07:32 2012 -0600 Work on factoring a little --- XS.xs | 32 ++++++++++++++++++++++++++------ factor.c | 56 +++++++++++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 81 insertions(+), 7 deletions(-) diff --git a/XS.xs b/XS.xs index 54937c0..fb0085b 100644 --- a/XS.xs +++ b/XS.xs @@ -227,26 +227,46 @@ _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 */ - /* For sufficiently large n, try more complex methods. */ /* SQUFOF (succeeds 98-99.9%) */ - split_success = squfof_factor(n, tofac_stack+ntofac, 256*1024)-1; - /* A few rounds of Pollard rho (succeeds most of the rest) */ if (!split_success) { - split_success = prho_factor(n, tofac_stack+ntofac, 400)-1; + 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) { 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); UV f = tlim; UV m = tlim % 30; UV limit = (UV) (sqrt(n)+0.1); @@ -317,7 +337,7 @@ fermat_factor(IN UV n, IN UV maxrounds = 64*1024*1024) SIMPLE_FACTOR(fermat_factor, n, maxrounds); void -holf_factor(IN UV n, IN UV maxrounds = 64*1024*1024) +holf_factor(IN UV n, IN UV maxrounds = 8*1024*1024) PPCODE: SIMPLE_FACTOR(holf_factor, n, maxrounds); @@ -337,7 +357,7 @@ 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 = 4*1024*1024) +pminus1_factor(IN UV n, IN UV maxrounds = 1*1024*1024) PPCODE: SIMPLE_FACTOR(pminus1_factor, n, maxrounds); diff --git a/factor.c b/factor.c index 42ad2d0..e8621cd 100644 --- a/factor.c +++ b/factor.c @@ -406,7 +406,7 @@ int holf_factor(UV n, UV *factors, UV rounds) MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in holf_factor"); for (i = 1; i <= rounds; i++) { - s = sqrt(n*i); /* TODO: overflow here */ + s = sqrt( (double)n * (double)i ); if ( (s*s) != (n*i) ) s++; m = sqrmod(s, n); if (is_perfect_square(m, &f)) { @@ -467,6 +467,7 @@ int pbrent_factor(UV n, UV *factors, UV rounds) * 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; @@ -502,6 +503,59 @@ int prho_factor(UV n, UV *factors, UV rounds) factors[0] = n; return 1; } +#else +int prho_factor(UV n, UV *factors, UV rounds) +{ + UV a, f, i, m, oldU, oldV; + const UV inner = 8; + UV U = 7; + UV V = 7; + + MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in prho_factor"); + + switch (n%8) { + case 1: a = 1; break; + case 3: a = 2; break; + case 5: a = 3; break; + case 7: a = 5; break; + default: a = 7; break; + } + + rounds = (rounds + inner - 1) / inner; + + while (rounds-- > 0) { + m = 1; oldU = U; oldV = V; + for (i = 0; i < inner; i++) { + U = sqraddmod(U, a, n); + V = sqraddmod(V, a, n); + V = sqraddmod(V, a, n); + f = (U > V) ? U-V : V-U; + m = mulmod(m, f, n); + } + f = gcd_ui(m, n); + if (f == 1) + continue; + if (f == n) { /* backup */ + U = oldU; V = oldV; + i = inner; + do { + 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); + } while (f == 1 && i-- != 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; +} +#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