Subject: factor: missed optimizations

From: Denys Vlasenko <[email-masked]>
Date: Sun, 20 Dec 2020 17:16:49 +0100

factor ---debug output is rather rudimentary.
In order to understand the logic, I added more printouts
(try attached debug patch)
and immediately noticed some obvious snafus:

$ src/factor ---debug 340282366920938461286658806734041124249
[using arbitrary-precision arithmetic] 340282366920938461286658806734041124249:[trial division 340282366920938461286658806734041124249] [trial stopped 4999(max)] [is number 34028236692093846128665880673404112424
9 prime?] [no] [pollard-rho (1)] 1:1:2:4:8:16:32:64:128:256:512:1024:2048:4096:8192:16384:32768:65536:131072:262144:524288:1048576:2097152:4194304:8388608:16777216:33554432:67108864:134217728:268435456:536870912
[factor found][rho factor 18446744073709551557] [trial division 18446744073709551556] [factor 2] [factor 2] [factor 11] [factor 137] [factor 547] [trial stopped 4999(max)] [is number 5594472617641 prime?] [trial
  division 5594472617640] [factor 2] [factor 2] [factor 2] [factor 3] [factor 5] [factor 1427] [factor 2131] [trial stopped 2131] [is number 15331 prime?] [yes] [factor 15331] [yes] ...

Look at the ens of the last line.
We don't need to check whether 15331 is prime. We know it is: we just stopped
trial division because next prime divisor to try is larger than
sqrt(number_we_try_to_factor) - therefore number_we_try_to_factor
has no more divisors, ergo it's prime. The fix:

    for (unsigned int i = 1; i <= PRIMES_PTAB_ENTRIES;)
      {
        if (! mpz_divisible_ui_p (t, p))
          {
            p += primes_diff[i++];
            if (mpz_cmp_ui (t, p * p) < 0)
+            mp_factor_insert (factors, t);
+            mpz_set_ui (t, 1);
              break;
          }
        else

... [factor 5594472617641] [factor
  18446744073709551557] [trial division 18446744073709551556] [factor 2] [factor 2] [factor 11] [factor 137] [factor 547] [trial stopped 4999(max)] [is number 5594472617641 prime?] [trial division 5594472617640]
[factor 2] [factor 2] [factor 2] [factor 3] [factor 5] [factor 1427] [factor 2131] [trial stopped 2131] [is number 15331 prime?] [yes] [factor 15331] [yes] [factor 5594472617641] [factor 18446744073709551557] [p
ollard-rho end]
   18446744073709551557 18446744073709551557

Second optimization is this: as you see, I deliberately fed it a square
(the factorization result is 18446744073709551557^2).
It happily did not notice it and spent an awful lot of time in pollard rho:

  1:1:2:4:8:16:32:64:128:256:512:1024:2048:4096:8192:16384:32768:65536:131072:262144:524288:1048576:2097152:4194304:8388608:16777216:33554432:67108864:134217728:268435456:536870912

A quick check for square drastically cuts down execution time:

  /* Use Pollard-rho to compute the prime factors of
     arbitrary-precision T, and put the results in FACTORS.  */
  static void
+mp_factor_big (mpz_t t, struct mp_factors *factors)
+{
+  /* the _big() function assumes trial division was already tried */
+  devmsg ("[is number prime?] ");
+  if (mp_prime_p (t))
+    mp_factor_insert (factors, t);
+  else
+    if (mpz_perfect_square_p (t))
+      {
+        struct mp_factors f2;
+        devmsg ("[it is square] ");
+        mpz_sqrt (t, t);
+        mp_factor_init (&f2);
+        mp_factor_big (t, &f2);
+        for (unsigned long int i = 0; i < f2.nfactors; i++)
+          {
+            unsigned long int e = f2.e[i];
+            while (e--) {                           //
+              mp_factor_insert (factors, f2.p[i]);  //FIXME: obvious optimization:
+              mp_factor_insert (factors, f2.p[i]);  // need to introduce/use mp_factor_insert_multiplicity(..., e*2);
+            }
+          }
+        mp_factor_clear (&f2);
+      }
+    else
+      mp_factor_using_pollard_rho (t, 1, factors);
+}
+
+static void
  mp_factor (mpz_t t, struct mp_factors *factors)
  {
    mp_factor_init (factors);
@@ -2265,11 +2309,7 @@ mp_factor (mpz_t t, struct mp_factors *f
.
        if (mpz_cmp_ui (t, 1) != 0)
          {
-          devmsg ("[is number prime?] ");
-          if (mp_prime_p (t))
-            mp_factor_insert (factors, t);
-          else
-            mp_factor_using_pollard_rho (t, 1, factors);
+          mp_factor_big (t, factors);
          }
      }
  }

diff -urp coreutils-8.32/src/factor.c coreutils-8.32-TT/src/factor.c
--- coreutils-8.32/src/factor.c	2020-01-01 15:13:12.000000000 +0100
+++ coreutils-8.32-TT/src/factor.c	2020-12-20 17:08:33.022081567 +0100
@@ -620,6 +620,8 @@ mp_factor_insert (struct mp_factors *fac
   unsigned long int *e  = factors->e;
   long i;
 
+gmp_printf ("[factor %Zd] ", prime);
+
   /* Locate position for insert new or increment e.  */
   for (i = nfactors - 1; i >= 0; i--)
     {
@@ -842,7 +844,8 @@ mp_factor_using_division (mpz_t t, struc
   mpz_t q;
   unsigned long int p;
 
-  devmsg ("[trial division] ");
+gmp_printf ("[trial division %Zd] ", t);
+//  devmsg ("[trial division] ");
 
   mpz_init (q);
 
@@ -861,6 +864,8 @@ mp_factor_using_division (mpz_t t, struc
         {
           p += primes_diff[i++];
           if (mpz_cmp_ui (t, p * p) < 0)
+            mp_factor_insert (factors, t);
+            mpz_set_ui (t, 1);
             break;
         }
       else
@@ -1407,6 +1412,7 @@ mp_prime_p (mpz_t n)
   if (!mp_millerrabin (n, nm1, a, tmp, q, k))
     {
       is_prime = false;
+gmp_printf ("[ML1:no] ");
       goto ret2;
     }
 
@@ -1414,6 +1420,7 @@ mp_prime_p (mpz_t n)
     {
       /* Factor n-1 for Lucas.  */
       mpz_set (tmp, nm1);
+gmp_printf ("[ML does not know, trying Lucas] ");
       mp_factor (tmp, &factors);
     }
 
@@ -1438,13 +1445,16 @@ mp_prime_p (mpz_t n)
         }
 
       if (is_prime)
+{gmp_printf ("[Lucas:yes] ");
         goto ret1;
+}
 
       mpz_add_ui (a, a, primes_diff[r]);        /* Establish new base.  */
 
       if (!mp_millerrabin (n, nm1, a, tmp, q, k))
         {
           is_prime = false;
+gmp_printf ("[ML2:no] ");
           goto ret1;
         }
     }
@@ -1692,6 +1702,7 @@ mp_factor_using_pollard_rho (mpz_t n, un
     {
       for (;;)
         {
+printf("%llu", k);
           do
             {
               mpz_mul (t, x, x);
@@ -1715,6 +1726,7 @@ mp_factor_using_pollard_rho (mpz_t n, un
           mpz_set (z, x);
           k = l;
           l = 2 * l;
+devmsg(":");
           for (unsigned long long int i = 0; i < k; i++)
             {
               mpz_mul (t, x, x);
@@ -1725,6 +1737,7 @@ mp_factor_using_pollard_rho (mpz_t n, un
         }
 
     factor_found:
+gmp_printf("[factor found]");
       do
         {
           mpz_mul (t, y, y);
@@ -1738,6 +1751,7 @@ mp_factor_using_pollard_rho (mpz_t n, un
 
       mpz_divexact (n, n, t);   /* divide by t, before t is overwritten */
 
+gmp_printf("[rho factor %Zd] ", t);
       if (!mp_prime_p (t))
         {
           devmsg ("[composite factor--restarting pollard-rho] ");
@@ -1759,6 +1773,7 @@ mp_factor_using_pollard_rho (mpz_t n, un
       mpz_mod (y, y, n);
     }
 
+  devmsg ("[pollard-rho end] ");
   mpz_clears (P, t2, t, z, x, y, NULL);
 }
 #endif
@@ -2255,6 +2270,35 @@ factor (uintmax_t t1, uintmax_t t0, stru
 /* Use Pollard-rho to compute the prime factors of
    arbitrary-precision T, and put the results in FACTORS.  */
 static void
+mp_factor_big (mpz_t t, struct mp_factors *factors)
+{
+  /* the _big() function assumes trial division was already tried */
+  devmsg ("[is number prime?] ");
+  if (mp_prime_p (t))
+    mp_factor_insert (factors, t);
+  else
+    if (mpz_perfect_square_p (t))
+      {
+        struct mp_factors f2;
+        devmsg ("[it is square] ");
+        mpz_sqrt (t, t);
+        mp_factor_init (&f2);
+        mp_factor_big (t, &f2);
+        for (unsigned long int i = 0; i < f2.nfactors; i++)
+          {
+            unsigned long int e = f2.e[i];
+            while (e--) {
+              mp_factor_insert (factors, f2.p[i]);
+              mp_factor_insert (factors, f2.p[i]);
+            }
+          }
+        mp_factor_clear (&f2);
+      }
+    else
+      mp_factor_using_pollard_rho (t, 1, factors);
+}
+
+static void
 mp_factor (mpz_t t, struct mp_factors *factors)
 {
   mp_factor_init (factors);
@@ -2265,11 +2309,7 @@ mp_factor (mpz_t t, struct mp_factors *f
 
       if (mpz_cmp_ui (t, 1) != 0)
         {
-          devmsg ("[is number prime?] ");
-          if (mp_prime_p (t))
-            mp_factor_insert (factors, t);
-          else
-            mp_factor_using_pollard_rho (t, 1, factors);
+          mp_factor_big (t, factors);
         }
     }
 }
@@ -2599,6 +2639,8 @@ do_stdin (void)
 int
 main (int argc, char **argv)
 {
+  setvbuf(stdout, NULL, _IONBF, 0);
+
   initialize_main (&argc, &argv);
   set_program_name (argv[0]);
   setlocale (LC_ALL, "");