[libmath-prime-util-perl] 47/59: Work on factoring a little

[Partha P. Mukherjee]

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 <dana at 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