Re: factor [Was: coreutils v5.2.1 - stat.c]

[ThMO]

Hello Eric, Paul and others listening,

attached you'll find my final patch for src/factor.c combining the two
different functions into one, as the underlying computation is the same,
except for the initial negation.

Additionally the error message regarding to factorize 0 will go to stderr,
but I do not return an error status in this case.

This combined function factorizes the two given examples within the
info-file somewhat faster than before (taking 96.4% of the previous time).

As far as the execution speed is concerned, it's IMHO not worth the trouble
for seeking a faster algorithm, since e.g. pari (v1.39) factorizes the above
two examples within microseconds.
FYI, even pari outputs -1 as the 1st factor for negative values.

Regarding to the naming of this command:
IMHO it should be named the same for all platforms as this makes much less
trouble when it comes to using `factor' inside shell-scripts, but then again
YMMV in this case.
Maybe it could be named factorize...?

THX for listening.

CU Tom.
(Thomas M.Ott)
Germany

--- coreutils-5.2.1/src/factor.c.orig   2004-01-21 23:27:02.000000000 +0100
+++ coreutils-5.2.1/src/factor.c        2005-10-02 14:02:25.000000000 +0200
@@ -43,6 +43,8 @@
 /* FIXME: if this program is ever modified to factor integers larger
    than 2^128, this constant (and the algorithm :-) will have to change.  */
 #define MAX_N_FACTORS 128
+#undef MAX_N_FACTORS
+#define        MAX_N_FACTORS   ( sizeof( uintmax_t)* CHAR_BIT- 1)
 
 /* The trial divisor increment wheel.  Use it to skip over divisors that
    are composites of 2, 3, 5, 7, or 11.  The part from WHEEL_START up to
@@ -91,78 +93,145 @@
   exit (status);
 }
 
-/* FIXME: comment */
+/* using a union simplifies cast orgies in order to omit warnings (ThMO) */
+
+union factor {
+    intmax_t   sfac;
+    uintmax_t  ufac;
+  };
+
+#if    HAVE_LONG_LONG
+  #define  ONE 1LL
+#else
+  #define  ONE 1L
+#endif
+
+/* sfactor(): factorize (un)signed long (long) ints -- (ThMO) */
 
 static int
-factor (uintmax_t n0, int max_n_factors, uintmax_t *factors)
+sfactor( const union factor sn0, const int is_unsigned, union factor *factors)
 {
-  register uintmax_t n = n0, d, q;
-  int n_factors = 0;
-  unsigned int const *w = wheel_tab;
-
-  if (n < 1)
-    return n_factors;
-
-  /* The exit condition in the following loop is correct because
-     any time it is tested one of these 3 conditions holds:
-      (1) d divides n
-      (2) n is prime
-      (3) n is composite but has no factors less than d.
-     If (1) or (2) obviously the right thing happens.
-     If (3), then since n is composite it is >= d^2. */
-
-  d = 2;
-  do
-    {
-      q = n / d;
-      while (n == q * d)
-       {
-         assert (n_factors < max_n_factors);
-         factors[n_factors++] = d;
-         n = q;
-         q = n / d;
-       }
-      d += *(w++);
-      if (w == WHEEL_END)
-       w = WHEEL_START;
-    }
-  while (d <= q);
-
-  if (n != 1 || n0 == 1)
-    {
-      assert (n_factors < max_n_factors);
-      factors[n_factors++] = n;
+  union factor        sn;
+  uintmax_t           n, d, q;
+  int                 n_factors;
+  const unsigned int  *w;
+
+  n_factors= 0;
+  sn= sn0;
+  if ( !is_unsigned && sn.sfac < 0)
+    {
+      factors[ n_factors++].sfac= -1;
+      if ( sn.sfac == ( ONE << ( sizeof( n)* CHAR_BIT- 1)))
+        {
+          factors[ n_factors++].ufac= 2;
+          sn.sfac >>= 1;
+        }
+      sn.sfac= -sn.sfac;
+    }
+
+  /* from now on the value is guaranteed to be positive, so force further
+   * calculations with unsigned quantities, since the evaluation of the
+   * remainder from a division will always have the correct sign, unlike
+   * the signed modulus operation - furthermore unsigned division is mostly
+   * somewhat faster than signed division -- (ThMO)
+   */
+  if ( sn.ufac <= 1)
+    return( n_factors);
+
+  for ( n= (uintmax_t) sn.sfac;  ( n & 1) == 0;  n >>= 1)
+    factors[ n_factors++].ufac= 2;     /* factorize by 2 the easy way... */
+
+  if ( n >= ( d= 3))
+    {
+      /* The exit condition in the following loop is correct because
+         any time it is tested one of these 3 conditions holds:
+          (1) d divides n
+          (2) n is prime
+          (3) n is composite but has no factors less than d.
+         If (1) or (2) obviously the right thing happens.
+         If (3), then since n is composite it is >= d^2. */
+
+      w= &wheel_tab[ 1];
+      do
+        {
+          while ( q= n/ d,  n == q* d)
+            {
+              factors[ n_factors++].ufac= d;
+              n= q;
+            }
+          d += *w++;
+          if ( w == WHEEL_END)
+            w -= WHEEL_END- WHEEL_START;
+        }
+      while ( d <= q);
+      if ( n > 1)
+        factors[ n_factors++].ufac= n;
     }
-
-  return n_factors;
+  return( n_factors);
 }
 
-/* FIXME: comment */
+/* print_factors(): try to factorize given number or state an
+ *                  appropriate error message
+ */
+
+#define        ary_size( t)    ( sizeof( (t))/ sizeof( (t)[ 0]) )
 
 static int
 print_factors (const char *s)
 {
-  uintmax_t factors[MAX_N_FACTORS];
-  uintmax_t n;
+  union factor  factors[ MAX_N_FACTORS+ 1], fn;
+  int           is_unsigned;
   int n_factors;
   int i;
-  char buf[INT_BUFSIZE_BOUND (uintmax_t)];
+  char buf[ 1+ INT_BUFSIZE_BOUND (uintmax_t)];
   strtol_error err;
 
-  if ((err = xstrtoumax (s, NULL, 10, &n, "")) != LONGINT_OK)
+  is_unsigned= 0;
+  if ( ( err= xstrtoimax( s, NULL, 10, &fn.sfac, "")) != LONGINT_OK)
     {
-      if (err == LONGINT_OVERFLOW)
-       error (0, 0, _("`%s' is too large"), s);
-      else
-       error (0, 0, _("`%s' is not a valid positive integer"), s);
-      return 1;
-    }
-  n_factors = factor (n, MAX_N_FACTORS, factors);
-  printf ("%s:", umaxtostr (n, buf));
-  for (i = 0; i < n_factors; i++)
-    printf (" %s", umaxtostr (factors[i], buf));
-  putchar ('\n');
-  return 0;
+      if ( err != LONGINT_OVERFLOW)
+        {
+          error( 0, 0, _( "%s is not a valid integer"), s);
+          return( 1);
+        }
+      else if ( fn.sfac < 0)   /* underflow? */
+        {
+          error( 0, 0, _( "%s is too small"), s);
+          return( 1);
+        }
+      /* overflowing a signed long (long) int -> try it unsigned then */
+      err= xstrtoumax( s, NULL, 10, &fn.ufac, "");
+      if ( err != LONGINT_OK)
+        {
+          error( 0, 0, err == LONGINT_OVERFLOW
+                         ? _( "%s it too large")
+                         : _( "%s is not a valid positive integer"), s);
+          return( 1);
+        }
+      ++is_unsigned;
+    }
+
+  if ( fn.sfac == 0)
+    {
+      fprintf( stderr, "0: %s\n", _( "cannot be factorized"));
+      return( 0);
+    }
+
+  /* factorize (un)signed long (long) ints */
+  n_factors= sfactor( fn, is_unsigned, factors);
+  assert( n_factors <= ary_size( factors));
+  printf( "%s:", !is_unsigned
+                   ? imaxtostr( fn.sfac, buf) : umaxtostr( fn.ufac, buf));
+  if ( n_factors > 0)
+    {
+      i= 0;
+      if ( !is_unsigned)       /* output maybe -1 */
+        printf( " %s", imaxtostr( factors[ i++].sfac, buf));
+      while ( i < n_factors)
+        printf( " %s", umaxtostr( factors[ i++].ufac, buf));
+    }
+  printf( "\n");
+  return( 0);
 }
 
 static int