[PATCH] factor: remove unreachable SQUFOF code at compile time

[Pádraig Brady]

It was a little confusing as to whether the SQUFOF algorithm was
enabled, and in fact there were no options available to enable it.
Therefore clarify the 3 configurable behaviors for the code to
3 defines at the top of the program, and only include the SQUFOF
code if enabled at compile time.

$ size src/factor-before
   text    data     bss
  93997    1412    2504
$ size src/factor-after
   text    data     bss
  87885    1404    2504

* src/factor.c: Only include the SQUFOF factor code
when enabled via the USE_SQUFOF define.
* doc/coreutils.texi (factor invocation): Update note about
factor limits, as we can factor 128 bit numbers without GMP.
---
 doc/coreutils.texi |  8 +++---
 src/factor.c       | 79 ++++++++++++++++++++++++++++--------------------------
 2 files changed, 45 insertions(+), 42 deletions(-)

diff --git a/doc/coreutils.texi b/doc/coreutils.texi
index c988aca..0359867 100644
--- a/doc/coreutils.texi
+++ b/doc/coreutils.texi
@@ -16855,7 +16855,7 @@ n=$(echo "$M8 * $M9" | bc)
 Similarly, factoring the eighth Fermat number @math{2^{256}+1} takes
 about 20 seconds on the same machine.
 
-Factoring large numbers is, in general, hard.  The Pollard Rho
+Factoring large numbers is, in general, hard.  The Pollard-Brent rho
 algorithm used by @command{factor} is particularly effective for
 numbers with relatively small factors.  If you wish to factor large
 numbers which do not have small factors (for example, numbers which
@@ -16863,9 +16863,9 @@ are the product of two large primes), other methods are 
far better.
 
 If @command{factor} is built without using GNU MP, only
 single-precision arithmetic is available, and so large numbers
-(typically @math{2^{64}} and above) will not be supported.  The 
single-precision
-code uses an algorithm which is designed for factoring smaller
-numbers.
+(typically @math{2^{128}} and above) will not be supported.
+The single-precision code uses an algorithm which is designed
+for factoring smaller numbers.
 
 @exitstatus
 
diff --git a/src/factor.c b/src/factor.c
index 1d7d7c8..e4ce4a5 100644
--- a/src/factor.c
+++ b/src/factor.c
@@ -50,8 +50,8 @@
     (2) Check the nature of any non-factored part using Miller-Rabin for
         detecting composites, and Lucas for detecting primes.
     (3) Factor any remaining composite part using the Pollard-Brent rho
-        algorithm or the SQUFOF algorithm, checking status of found factors
-        again using Miller-Rabin and Lucas.
+        algorithm or if USE_SQUFOF is defined to 1, try that first.
+        Status of found factors are checked again using Miller-Rabin and Lucas.
 
     We prefer using Hensel norm in the divisions, not the more familiar
     Euclidian norm, since the former leads to much faster code.  In the
@@ -84,6 +84,23 @@
       the -w option.
 */
 
+/* Whether to recursively factor to prove primality,
+   or run faster probabilistic tests.  */
+#ifndef PROVE_PRIMALITY
+# define PROVE_PRIMALITY 1
+#endif
+
+/* Faster for certain ranges but less general.  */
+#ifndef USE_SQUFOF
+# define USE_SQUFOF 0
+#endif
+
+/* Output SQUFOF statistics.  */
+#ifndef STAT_SQUFOF
+# define STAT_SQUFOF 0
+#endif
+
+
 #include <config.h>
 #include <getopt.h>
 #include <stdio.h>
@@ -114,10 +131,6 @@
 /* Token delimiters when reading from a file.  */
 #define DELIM "\n\t "
 
-#ifndef STAT_SQUFOF
-# define STAT_SQUFOF 0
-#endif
-
 #ifndef USE_LONGLONG_H
 /* With the way we use longlong.h, it's only safe to use
    when UWtype = UHWtype, as there were various cases
@@ -208,10 +221,6 @@ const unsigned char factor_clz_tab[129] =
    FIXME: this is just an ugly band-aid.  Fix it properly.  */
 #endif
 
-enum alg_type { ALG_POLLARD_RHO = 1, ALG_SQUFOF = 2 };
-
-static enum alg_type alg;
-
 /* 2*3*5*7*11...*101 is 128 bits, and has 26 prime factors */
 #define MAX_NFACTS 26
 
@@ -691,7 +700,7 @@ verify (W <= WIDE_UINT_BITS);
 static bool dev_debug = false;
 
 /* Prove primality or run probabilistic tests.  */
-static bool flag_prove_primality = true;
+static bool flag_prove_primality = PROVE_PRIMALITY;
 
 /* Number of Miller-Rabin tests to run when not proving primality. */
 #define MR_REPS 25
@@ -1741,6 +1750,7 @@ mp_factor_using_pollard_rho (mpz_t n, unsigned long int a,
 }
 #endif
 
+#if USE_SQUFOF
 /* FIXME: Maybe better to use an iteration converging to 1/sqrt(n)?  If
    algorithm is replaced, consider also returning the remainder. */
 static uintmax_t _GL_ATTRIBUTE_CONST
@@ -1812,10 +1822,10 @@ isqrt2 (uintmax_t nh, uintmax_t nl)
 }
 
 /* MAGIC[N] has a bit i set iff i is a quadratic residue mod N. */
-#define MAGIC64 0x0202021202030213ULL
-#define MAGIC63 0x0402483012450293ULL
-#define MAGIC65 0x218a019866014613ULL
-#define MAGIC11 0x23b
+# define MAGIC64 0x0202021202030213ULL
+# define MAGIC63 0x0402483012450293ULL
+# define MAGIC65 0x218a019866014613ULL
+# define MAGIC11 0x23b
 
 /* Return the square root if the input is a square, otherwise 0. */
 static uintmax_t _GL_ATTRIBUTE_CONST
@@ -1860,7 +1870,7 @@ static const unsigned short invtab[0x81] =
 
 /* Compute q = [u/d], r = u mod d.  Avoids slow hardware division for the case
    that q < 0x40; here it instead uses a table of (Euclidian) inverses.  */
-#define div_smallq(q, r, u, d)                                          \
+# define div_smallq(q, r, u, d)                                          \
   do {                                                                  \
     if ((u) / 0x40 < (d))                                               \
       {                                                                 \
@@ -1928,7 +1938,8 @@ static const unsigned short invtab[0x81] =
    0.9295  7        =    7 = 3
    0.9934  11       =   11 = 3
 */
-#define QUEUE_SIZE 50
+# define QUEUE_SIZE 50
+#endif
 
 #if STAT_SQUFOF
 # define Q_FREQ_SIZE 50
@@ -1937,6 +1948,7 @@ static unsigned int q_freq[Q_FREQ_SIZE + 1];
 # define MIN(a,b) ((a) < (b) ? (a) : (b))
 #endif
 
+#if USE_SQUFOF
 /* Return true on success.  Expected to fail only for numbers
    >= 2^{2*W_TYPE_SIZE - 2}, or close to that limit. */
 static bool
@@ -2058,10 +2070,10 @@ factor_using_squfof (uintmax_t n1, uintmax_t n0, struct 
factors *factors)
 
           IF_LINT (assert (q > 0 && Q > 0));
 
-#if STAT_SQUFOF
+# if STAT_SQUFOF
           q_freq[0]++;
           q_freq[MIN (q, Q_FREQ_SIZE)]++;
-#endif
+# endif
 
           if (Q <= L1)
             {
@@ -2145,10 +2157,10 @@ factor_using_squfof (uintmax_t n1, uintmax_t n0, struct 
factors *factors)
                       div_smallq (q, rem, S+P, Q);
                       P1 = S - rem;     /* P1 = q*Q - P */
 
-#if STAT_SQUFOF
+# if STAT_SQUFOF
                       q_freq[0]++;
                       q_freq[MIN (q, Q_FREQ_SIZE)]++;
-#endif
+# endif
                       if (P == P1)
                         break;
                       t = Q1 + q * (P - P1);
@@ -2192,9 +2204,10 @@ factor_using_squfof (uintmax_t n1, uintmax_t n0, struct 
factors *factors)
     }
   return false;
 }
+#endif
 
 /* Compute the prime factors of the 128-bit number (T1,T0), and put the
-   results in FACTORS.  Use the algorithm selected by the global ALG.  */
+   results in FACTORS.  */
 static void
 factor (uintmax_t t1, uintmax_t t0, struct factors *factors)
 {
@@ -2213,9 +2226,10 @@ factor (uintmax_t t1, uintmax_t t0, struct factors 
*factors)
     factor_insert_large (factors, t1, t0);
   else
     {
-      if (alg == ALG_SQUFOF)
-        if (factor_using_squfof (t1, t0, factors))
-          return;
+#if USE_SQUFOF
+      if (factor_using_squfof (t1, t0, factors))
+        return;
+#endif
 
       if (t1 == 0)
         factor_using_pollard_rho (t0, 1, factors);
@@ -2576,8 +2590,6 @@ main (int argc, char **argv)
   atexit (close_stdout);
   atexit (lbuf_flush);
 
-  alg = ALG_POLLARD_RHO;        /* Default to Pollard rho */
-
   int c;
   while ((c = getopt_long (argc, argv, "", long_options, NULL)) != -1)
     {
@@ -2587,14 +2599,6 @@ main (int argc, char **argv)
           dev_debug = true;
           break;
 
-        case 's':
-          alg = ALG_SQUFOF;
-          break;
-
-        case 'w':
-          flag_prove_primality = false;
-          break;
-
         case_GETOPT_HELP_CHAR;
 
         case_GETOPT_VERSION_CHAR (PROGRAM_NAME, AUTHORS);
@@ -2605,8 +2609,7 @@ main (int argc, char **argv)
     }
 
 #if STAT_SQUFOF
-  if (alg == ALG_SQUFOF)
-    memset (q_freq, 0, sizeof (q_freq));
+  memset (q_freq, 0, sizeof (q_freq));
 #endif
 
   bool ok;
@@ -2621,7 +2624,7 @@ main (int argc, char **argv)
     }
 
 #if STAT_SQUFOF
-  if (alg == ALG_SQUFOF && q_freq[0] > 0)
+  if (q_freq[0] > 0)
     {
       double acc_f;
       printf ("q  freq.  cum. freq.(total: %d)\n", q_freq[0]);
-- 
2.5.0