Merge branch 'busybox' into merge

1 files changed, 132 insertions, 77 deletions

diff --git a/coreutils/factor.c b/coreutils/factor.c
index 47fe179dc..a7a5a5030 100644
--- a/coreutils/factor.c
+++ b/coreutils/factor.c
@@ -19,6 +19,7 @@
 //usage:       "Print prime factors"
 
 #include "libbb.h"
+#include "common_bufsiz.h"
 
 #if 0
 # define dbg(...) bb_error_msg(__VA_ARGS__)
@@ -42,9 +43,116 @@ typedef unsigned long half_t;
 #error Cant find an integer type which is half as wide as ullong
 #endif
 
+/* The trial divisor increment wheel.  Use it to skip over divisors that
+ * are composites of 2, 3, 5, 7, or 11.
+ * Larger wheels improve sieving only slightly, but quickly grow in size
+ * (adding just one prime, 13, results in 5766 element sieve).
+ */
+#define R(a,b,c,d,e,f,g,h,i,j,A,B,C,D,E,F,G,H,I,J) \
+        (((uint64_t)(a<<0) | (b<<3) | (c<<6) | (d<<9) | (e<<12) | (f<<15) | (g<<18) | (h<<21) | (i<<24) | (j<<27)) << 1) | \
+        (((uint64_t)(A<<0) | (B<<3) | (C<<6) | (D<<9) | (E<<12) | (F<<15) | (G<<18) | (H<<21) | (I<<24) | (J<<27)) << 31)
+#define P(a,b,c,d,e,f,g,h,i,j,A,B,C,D,E,F,G,H,I,J) \
+        R(      (a/2),(b/2),(c/2),(d/2),(e/2),(f/2),(g/2),(h/2),(i/2),(j/2), \
+                (A/2),(B/2),(C/2),(D/2),(E/2),(F/2),(G/2),(H/2),(I/2),(J/2)  )
+static const uint64_t packed_wheel[] = {
+        /*1, 2, 2, 4, 2,*/
+        P( 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4), //01
+        P( 2, 4, 2, 4,14, 4, 6, 2,10, 2, 6, 6, 4, 2, 4, 6, 2,10, 2, 4), //02
+        P( 2,12,10, 2, 4, 2, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4), //03
+        P( 6, 8, 4, 2, 4, 6, 8, 6,10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2), //04
+        P( 6, 4, 2, 6,10, 2,10, 2, 4, 2, 4, 6, 8, 4, 2, 4,12, 2, 6, 4), //05
+        P( 2, 6, 4, 6,12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6,10, 2), //06
+        P( 4, 6, 2, 6, 4, 2, 4, 2,10, 2,10, 2, 4, 6, 6, 2, 6, 6, 4, 6), //07
+        P( 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6), //08
+        P( 8, 6, 4, 2,10, 2, 6, 4, 2, 4, 2,10, 2,10, 2, 4, 2, 4, 8, 6), //09
+        P( 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4), //10
+        P( 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2,10, 2), //11
+        P( 6, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6,10, 8, 4, 2, 4, 2), //12
+        P( 4, 8,10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 6, 2,10, 2), //13
+        P(10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8), //14
+        P( 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4), //15
+        P( 2, 4, 2,10, 2,10, 2, 4, 2, 4, 6, 2,10, 2, 4, 6, 8, 6, 4, 2), //16
+        P( 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 6), //17
+        P( 6, 2, 6, 6, 4, 2,10, 2,10, 2, 4, 2, 4, 6, 2, 6, 4, 2,10, 6), //18
+        P( 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2,12, 6, 4, 6, 2, 4, 6, 2), //19
+        P(12, 4, 2, 4, 8, 6, 4, 2, 4, 2,10, 2,10, 6, 2, 4, 6, 2, 6, 4), //20
+        P( 2, 4, 6, 6, 2, 6, 4, 2,10, 6, 8, 6, 4, 2, 4, 8, 6, 4, 6, 2), //21
+        P( 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 2, 4, 2,10,12, 2, 4, 2,10), //22
+        P( 2, 6, 4, 2, 4, 6, 6, 2,10, 2, 6, 4,14, 4, 2, 4, 2, 4, 8, 6), //23
+        P( 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4,12, 2,12), //24
+};
+#undef P
+#undef R
+#define WHEEL_START 5
+#define WHEEL_SIZE (5 + 24 * 20)
+#define square_count (((uint8_t*)&bb_common_bufsiz1)[0])
+#define wheel_tab    (((uint8_t*)&bb_common_bufsiz1) + 1)
+/*
+ * Why, you ask?
+ * plain byte array:
+ *      function                old     new   delta
+ *      wheel_tab                 -     485    +485
+ * 3-bit-packed insanity:
+ *      packed_wheel              -     192    +192
+ *      factor_main             108     171     +63
+ */
+static void unpack_wheel(void)
+{
+        int i;
+        uint8_t *p;
+
+        setup_common_bufsiz();
+        wheel_tab[0] = 1;
+        wheel_tab[1] = 2;
+        wheel_tab[2] = 2;
+        wheel_tab[3] = 4;
+        wheel_tab[4] = 2;
+        p = &wheel_tab[5];
+        for (i = 0; i < ARRAY_SIZE(packed_wheel); i++) {
+                uint64_t v = packed_wheel[i];
+                while ((v & 0xe) != 0) {
+                        *p = v & 0xe;
+                        //printf("%2u,", *p);
+                        p++;
+                        v >>= 3;
+                }
+                //printf("\n");
+        }
+}
+
+/* Prevent inlining, factorize() needs all help it can get with reducing register pressure */
+static NOINLINE void print_w(wide_t n)
+{
+        unsigned rep = square_count;
+        do
+                printf(" %"LL_FMT"u", n);
+        while (--rep != 0);
+}
+static NOINLINE void print_h(half_t n)
+{
+        print_w(n);
+}
+
+static void factorize(wide_t N);
+
 static half_t isqrt_odd(wide_t N)
 {
         half_t s = isqrt(N);
+        /* s^2 is <= N, (s+1)^2 > N */
+
+        /* If s^2 in fact is EQUAL to N, it's very lucky.
+         * Examples:
+         * factor 18446743988964486098 = 2 * 3037000493 * 3037000493
+         * factor 18446743902517389507 = 3 * 2479700513 * 2479700513
+         */
+        if ((wide_t)s * s == N) {
+                /* factorize sqrt(N), printing each factor twice */
+                square_count *= 2;
+                factorize(s);
+                /* Let caller know we recursed */
+                return 0;
+        }
+
         /* Subtract 1 from even s, odd s won't change: */
         /* (doesnt work for zero, but we know that s != 0 here) */
         s = (s - 1) | 1;
@@ -53,43 +161,20 @@ static half_t isqrt_odd(wide_t N)
 
 static NOINLINE void factorize(wide_t N)
 {
+        unsigned w;
         half_t factor;
         half_t max_factor;
-        // unsigned count3;
-        // unsigned count5;
-        // unsigned count7;
-        // ^^^^^^^^^^^^^^^ commented-out simple sieving code (easier to grasp).
-        // Faster sieving, using one word for potentially up to 6 counters:
-        // count upwards in each mask, counter "triggers" when it sets its mask to "100[0]..."
-        // 10987654321098765432109876543210 - bits 31-0 in 32-bit word
-        //    17777713333311111777775555333 - bit masks for counters for primes 3,5,7,11,13,17
-        //         100000100001000010001001 - value for adding 1 to each mask
-        //    10000010000010000100001000100 - value for checking that any mask reached msb
-        enum {
-                SHIFT_3 = 1 << 0,
-                SHIFT_5 = 1 << 3,
-                SHIFT_7 = 1 << 7,
-                INCREMENT_EACH = SHIFT_3 | SHIFT_5 | SHIFT_7,
-                MULTIPLE_OF_3 = 1 << 2,
-                MULTIPLE_OF_5 = 1 << 6,
-                MULTIPLE_OF_7 = 1 << 11,
-                MULTIPLE_DETECTED = MULTIPLE_OF_3 | MULTIPLE_OF_5 | MULTIPLE_OF_7,
-        };
-        unsigned sieve_word;
 
         if (N < 4)
                 goto end;
 
-        while (!(N & 1)) {
-                printf(" 2");
-                N >>= 1;
-        }
-
         /* The code needs to be optimized for the case where
          * there are large prime factors. For example,
          * this is not hard:
          * 8262075252869367027 = 3 7 17 23 47 101 113 127 131 137 823
-         * (the largest factor to test is only ~sqrt(823) = 28)
+         *  (the largest divisor to test for largest factor 823
+         *  is only ~sqrt(823) = 28, the entire factorization needs
+         *  only ~33 trial divisions)
          * but this is:
          * 18446744073709551601 = 53 348051774975651917
          * the last factor requires testing up to
@@ -98,70 +183,37 @@ static NOINLINE void factorize(wide_t N)
          * factor 18446744073709551557 (0xffffffffffffffc5).
          */
         max_factor = isqrt_odd(N);
-        // count3 = 3;
-        // count5 = 6;
-        // count7 = 9;
-        sieve_word = 0
-                /* initial count for SHIFT_n is (n-1)/2*3: */
-                + (MULTIPLE_OF_3 - 3 * SHIFT_3)
-                + (MULTIPLE_OF_5 - 6 * SHIFT_5)
-                + (MULTIPLE_OF_7 - 9 * SHIFT_7)
-                //+ (MULTIPLE_OF_11 - 15 * SHIFT_11)
-                //+ (MULTIPLE_OF_13 - 18 * SHIFT_13)
-                //+ (MULTIPLE_OF_17 - 24 * SHIFT_17)
-        ;
-        factor = 3;
+        if (!max_factor)
+                return; /* square was detected and recursively factored */
+        factor = 2;
+        w = 0;
         for (;;) {
+                half_t fw;
+
                 /* The division is the most costly part of the loop.
                  * On 64bit CPUs, takes at best 12 cycles, often ~20.
                  */
                 while ((N % factor) == 0) { /* not likely */
                         N = N / factor;
-                        printf(" %"HALF_FMT"u", factor);
+                        print_h(factor);
                         max_factor = isqrt_odd(N);
+                        if (!max_factor)
+                                return; /* square was detected */
                 }
- next_factor:
                 if (factor >= max_factor)
                         break;
-                factor += 2;
-                /* Rudimentary wheel sieving: skip multiples of 3, 5 and 7:
-                 * Every third odd number is divisible by three and thus isn't a prime:
-                 * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47...
-                 * ^ ^   ^  ^     ^  ^     ^  _     ^  ^     _  ^     ^  ^     ^
-                 * (^ = primes, _ = would-be-primes-if-not-divisible-by-5)
-                 * The numbers with space under them are excluded by sieve 3.
-                 */
-                // count7--;
-                // count5--;
-                // count3--;
-                // if (count3 && count5 && count7)
-                //      continue;
-                sieve_word += INCREMENT_EACH;
-                if (!(sieve_word & MULTIPLE_DETECTED))
+                fw = factor + wheel_tab[w];
+                if (fw < factor)
+                        break; /* overflow */
+                factor = fw;
+                w++;
+                if (w < WHEEL_SIZE)
                         continue;
-                /*
-                 * "factor" is multiple of 3 33% of the time (count3 reached 0),
-                 * else, multiple of 5 13% of the time,
-                 * else, multiple of 7 7.6% of the time.
-                 * Cumulatively, with 3,5,7 sieving we are here 54.3% of the time.
-                 */
-                // if (count3 == 0)
-                //      count3 = 3;
-                if (sieve_word & MULTIPLE_OF_3)
-                        sieve_word -= SHIFT_3 * 3;
-                // if (count5 == 0)
-                //      count5 = 5;
-                if (sieve_word & MULTIPLE_OF_5)
-                        sieve_word -= SHIFT_5 * 5;
-                // if (count7 == 0)
-                //      count7 = 7;
-                if (sieve_word & MULTIPLE_OF_7)
-                        sieve_word -= SHIFT_7 * 7;
-                goto next_factor;
+                w = WHEEL_START;
         }
  end:
         if (N > 1)
-                printf(" %"LL_FMT"u", N);
+                print_w(N);
         bb_putchar('\n');
 }
 
@@ -176,12 +228,15 @@ static void factorize_numstr(const char *numstr)
         if (errno)
                 bb_show_usage();
         printf("%"LL_FMT"u:", N);
+        square_count = 1;
         factorize(N);
 }
 
 int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
 int factor_main(int argc UNUSED_PARAM, char **argv)
 {
+        unpack_wheel();
+
         //// coreutils has undocumented option ---debug (three dashes)
         //getopt32(argv, "");
         //argv += optind;