ntpd: default to FEATURE_NTP_AUTH=y
[oweals/busybox.git] / coreutils / factor.c
index 0f838a735da9340d754e379ff4106e0133ca5f16..467e23a4d81aebbda0fcba7cab6f31d05f14c0d8 100644 (file)
  * Licensed under GPLv2, see file LICENSE in this source tree.
  */
 //config:config FACTOR
-//config:      bool "factor"
+//config:      bool "factor (2.6 kb)"
 //config:      default y
 //config:      help
-//config:        factor factorizes integers
+//config:      factor factorizes integers
 
 //applet:IF_FACTOR(APPLET(factor, BB_DIR_USR_BIN, BB_SUID_DROP))
 
 //kbuild:lib-$(CONFIG_FACTOR) += factor.o
 
 //usage:#define factor_trivial_usage
-//usage:       "NUMBER..."
+//usage:       "[NUMBER]..."
 //usage:#define factor_full_usage "\n\n"
 //usage:       "Print prime factors"
 
 #include "libbb.h"
 
+#if 0
+# define dbg(...) bb_error_msg(__VA_ARGS__)
+#else
+# define dbg(...) ((void)0)
+#endif
+
+typedef unsigned long long wide_t;
+
 #if ULLONG_MAX == (UINT_MAX * UINT_MAX + 2 * UINT_MAX)
 /* "unsigned" is half as wide as ullong */
 typedef unsigned half_t;
+#define HALF_MAX UINT_MAX
 #define HALF_FMT ""
 #elif ULLONG_MAX == (ULONG_MAX * ULONG_MAX + 2 * ULONG_MAX)
 /* long is half as wide as ullong */
 typedef unsigned long half_t;
+#define HALF_MAX ULONG_MAX
 #define HALF_FMT "l"
 #else
 #error Cant find an integer type which is half as wide as ullong
 #endif
 
-static void factorize(const char *numstr)
+static half_t isqrt_odd(wide_t N)
 {
-       unsigned long long N, factor2;
-       half_t factor;
-       unsigned count3;
-
-       /* Coreutils compat */
-       numstr = skip_whitespace(numstr);
-       if (*numstr == '+')
-               numstr++;
-
-       N = bb_strtoull(numstr, NULL, 10);
-       if (errno)
-               bb_show_usage();
+       half_t s = isqrt(N);
+       /* 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;
+       return s;
+}
 
-       printf("%llu:", N);
+static NOINLINE void factorize(wide_t N)
+{
+       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;
        }
 
-       count3 = 3;
+       /* 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)
+        * but this is:
+        * 18446744073709551601 = 53 348051774975651917
+        * the last factor requires testing up to
+        * 589959129 - about 100 million iterations.
+        * The slowest case (largest prime) for N < 2^64 is
+        * 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;
-       factor2 = 3 * 3;
        for (;;) {
-               unsigned long long nfactor2;
-
-               while ((N % factor) == 0) {
+               /* 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(" %u"HALF_FMT"", factor);
+                       printf(" %"HALF_FMT"u", factor);
+                       max_factor = isqrt_odd(N);
                }
  next_factor:
-               /* (f + 2)^2 = f^2 + 4*f + 4 = f^2 + 4*(f+1) */
-               nfactor2 = factor2 + 4 * (factor + 1);
-               if (nfactor2 < factor2) /* overflow? */
-                       break;
-               factor2 = nfactor2;
-               if (factor2 > N)
+               if (factor >= max_factor)
                        break;
                factor += 2;
-               /* Rudimentary wheel sieving: skip multiples of 3:
+               /* 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...
-                * ^  ^     ^   ^       ^  ^     ^  _     ^  ^     _  ^ (^primes)
+                * 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.
                 */
-               --count3;
-               if (count3 == 0) {
-                       count3 = 3;
-                       goto next_factor;
-               }
+               // count7--;
+               // count5--;
+               // count3--;
+               // if (count3 && count5 && count7)
+               //      continue;
+               sieve_word += INCREMENT_EACH;
+               if (!(sieve_word & MULTIPLE_DETECTED))
+                       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;
        }
  end:
        if (N > 1)
@@ -92,6 +165,20 @@ static void factorize(const char *numstr)
        bb_putchar('\n');
 }
 
+static void factorize_numstr(const char *numstr)
+{
+       wide_t N;
+
+       /* Leading + is ok (coreutils compat) */
+       if (*numstr == '+')
+               numstr++;
+       N = bb_strtoull(numstr, NULL, 10);
+       if (errno)
+               bb_show_usage();
+       printf("%llu:", N);
+       factorize(N);
+}
+
 int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
 int factor_main(int argc UNUSED_PARAM, char **argv)
 {
@@ -100,13 +187,33 @@ int factor_main(int argc UNUSED_PARAM, char **argv)
        //argv += optind;
        argv++;
 
-       if (!*argv)
-               //TODO: read from stdin
-               bb_show_usage();
+       if (!*argv) {
+               /* Read from stdin, several numbers per line are accepted */
+               for (;;) {
+                       char *numstr, *line;
+                       line = xmalloc_fgetline(stdin);
+                       if (!line)
+                               return EXIT_SUCCESS;
+                       numstr = line;
+                       for (;;) {
+                               char *end;
+                               numstr = skip_whitespace(numstr);
+                               if (!numstr[0])
+                                       break;
+                               end = skip_non_whitespace(numstr);
+                               if (*end != '\0')
+                                       *end++ = '\0';
+                               factorize_numstr(numstr);
+                               numstr = end;
+                       }
+                       free(line);
+               }
+       }
 
        do {
-               factorize(*argv);
+               /* Leading spaces are ok (coreutils compat) */
+               factorize_numstr(skip_whitespace(*argv));
        } while (*++argv);
 
-       fflush_stdout_and_exit(EXIT_SUCCESS);
+       return EXIT_SUCCESS;
 }