coreutils/factor.c

   1 /*
   2  * Copyright (C) 2017 Denys Vlasenko <vda.linux@googlemail.com>
   3  *
   4  * Licensed under GPLv2, see file LICENSE in this source tree.
   5  */
   6 //config:config FACTOR
   7 //config:       bool "factor (2.6 kb)"
   8 //config:       default y
   9 //config:       help
  10 //config:       factor factorizes integers
  11
  12 //applet:IF_FACTOR(APPLET(factor, BB_DIR_USR_BIN, BB_SUID_DROP))
  13
  14 //kbuild:lib-$(CONFIG_FACTOR) += factor.o
  15
  16 //usage:#define factor_trivial_usage
  17 //usage:       "[NUMBER]..."
  18 //usage:#define factor_full_usage "\n\n"
  19 //usage:       "Print prime factors"
  20
  21 #include "libbb.h"
  22
  23 #if 0
  24 # define dbg(...) bb_error_msg(__VA_ARGS__)
  25 #else
  26 # define dbg(...) ((void)0)
  27 #endif
  28
  29 typedef unsigned long long wide_t;
  30
  31 #if ULLONG_MAX == (UINT_MAX * UINT_MAX + 2 * UINT_MAX)
  32 /* "unsigned" is half as wide as ullong */
  33 typedef unsigned half_t;
  34 #define HALF_MAX UINT_MAX
  35 #define HALF_FMT ""
  36 #elif ULLONG_MAX == (ULONG_MAX * ULONG_MAX + 2 * ULONG_MAX)
  37 /* long is half as wide as ullong */
  38 typedef unsigned long half_t;
  39 #define HALF_MAX ULONG_MAX
  40 #define HALF_FMT "l"
  41 #else
  42 #error Cant find an integer type which is half as wide as ullong
  43 #endif
  44
  45 static half_t isqrt_odd(wide_t N)
  46 {
  47         half_t s = isqrt(N);
  48         /* Subtract 1 from even s, odd s won't change: */
  49         /* (doesnt work for zero, but we know that s != 0 here) */
  50         s = (s - 1) | 1;
  51         return s;
  52 }
  53
  54 static NOINLINE void factorize(wide_t N)
  55 {
  56         half_t factor;
  57         half_t max_factor;
  58         // unsigned count3;
  59         // unsigned count5;
  60         // unsigned count7;
  61         // ^^^^^^^^^^^^^^^ commented-out simple sieving code (easier to grasp).
  62         // Faster sieving, using one word for potentially up to 6 counters:
  63         // count upwards in each mask, counter "triggers" when it sets its mask to "100[0]..."
  64         // 10987654321098765432109876543210 - bits 31-0 in 32-bit word
  65         //    17777713333311111777775555333 - bit masks for counters for primes 3,5,7,11,13,17
  66         //         100000100001000010001001 - value for adding 1 to each mask
  67         //    10000010000010000100001000100 - value for checking that any mask reached msb
  68         enum {
  69                 SHIFT_3 = 1 << 0,
  70                 SHIFT_5 = 1 << 3,
  71                 SHIFT_7 = 1 << 7,
  72                 INCREMENT_EACH = SHIFT_3 | SHIFT_5 | SHIFT_7,
  73                 MULTIPLE_OF_3 = 1 << 2,
  74                 MULTIPLE_OF_5 = 1 << 6,
  75                 MULTIPLE_OF_7 = 1 << 11,
  76                 MULTIPLE_DETECTED = MULTIPLE_OF_3 | MULTIPLE_OF_5 | MULTIPLE_OF_7,
  77         };
  78         unsigned sieve_word;
  79
  80         if (N < 4)
  81                 goto end;
  82
  83         while (!(N & 1)) {
  84                 printf(" 2");
  85                 N >>= 1;
  86         }
  87
  88         /* The code needs to be optimized for the case where
  89          * there are large prime factors. For example,
  90          * this is not hard:
  91          * 8262075252869367027 = 3 7 17 23 47 101 113 127 131 137 823
  92          * (the largest factor to test is only ~sqrt(823) = 28)
  93          * but this is:
  94          * 18446744073709551601 = 53 348051774975651917
  95          * the last factor requires testing up to
  96          * 589959129 - about 100 million iterations.
  97          * The slowest case (largest prime) for N < 2^64 is
  98          * factor 18446744073709551557 (0xffffffffffffffc5).
  99          */
 100         max_factor = isqrt_odd(N);
 101         // count3 = 3;
 102         // count5 = 6;
 103         // count7 = 9;
 104         sieve_word = 0
 105                 /* initial count for SHIFT_n is (n-1)/2*3: */
 106                 + (MULTIPLE_OF_3 - 3 * SHIFT_3)
 107                 + (MULTIPLE_OF_5 - 6 * SHIFT_5)
 108                 + (MULTIPLE_OF_7 - 9 * SHIFT_7)
 109                 //+ (MULTIPLE_OF_11 - 15 * SHIFT_11)
 110                 //+ (MULTIPLE_OF_13 - 18 * SHIFT_13)
 111                 //+ (MULTIPLE_OF_17 - 24 * SHIFT_17)
 112         ;
 113         factor = 3;
 114         for (;;) {
 115                 /* The division is the most costly part of the loop.
 116                  * On 64bit CPUs, takes at best 12 cycles, often ~20.
 117                  */
 118                 while ((N % factor) == 0) { /* not likely */
 119                         N = N / factor;
 120                         printf(" %"HALF_FMT"u", factor);
 121                         max_factor = isqrt_odd(N);
 122                 }
 123  next_factor:
 124                 if (factor >= max_factor)
 125                         break;
 126                 factor += 2;
 127                 /* Rudimentary wheel sieving: skip multiples of 3, 5 and 7:
 128                  * Every third odd number is divisible by three and thus isn't a prime:
 129                  * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47...
 130                  * ^ ^   ^  ^     ^  ^     ^  _     ^  ^     _  ^     ^  ^     ^
 131                  * (^ = primes, _ = would-be-primes-if-not-divisible-by-5)
 132                  * The numbers with space under them are excluded by sieve 3.
 133                  */
 134                 // count7--;
 135                 // count5--;
 136                 // count3--;
 137                 // if (count3 && count5 && count7)
 138                 //      continue;
 139                 sieve_word += INCREMENT_EACH;
 140                 if (!(sieve_word & MULTIPLE_DETECTED))
 141                         continue;
 142                 /*
 143                  * "factor" is multiple of 3 33% of the time (count3 reached 0),
 144                  * else, multiple of 5 13% of the time,
 145                  * else, multiple of 7 7.6% of the time.
 146                  * Cumulatively, with 3,5,7 sieving we are here 54.3% of the time.
 147                  */
 148                 // if (count3 == 0)
 149                 //      count3 = 3;
 150                 if (sieve_word & MULTIPLE_OF_3)
 151                         sieve_word -= SHIFT_3 * 3;
 152                 // if (count5 == 0)
 153                 //      count5 = 5;
 154                 if (sieve_word & MULTIPLE_OF_5)
 155                         sieve_word -= SHIFT_5 * 5;
 156                 // if (count7 == 0)
 157                 //      count7 = 7;
 158                 if (sieve_word & MULTIPLE_OF_7)
 159                         sieve_word -= SHIFT_7 * 7;
 160                 goto next_factor;
 161         }
 162  end:
 163         if (N > 1)
 164                 printf(" %llu", N);
 165         bb_putchar('\n');
 166 }
 167
 168 static void factorize_numstr(const char *numstr)
 169 {
 170         wide_t N;
 171
 172         /* Leading + is ok (coreutils compat) */
 173         if (*numstr == '+')
 174                 numstr++;
 175         N = bb_strtoull(numstr, NULL, 10);
 176         if (errno)
 177                 bb_show_usage();
 178         printf("%llu:", N);
 179         factorize(N);
 180 }
 181
 182 int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
 183 int factor_main(int argc UNUSED_PARAM, char **argv)
 184 {
 185         //// coreutils has undocumented option ---debug (three dashes)
 186         //getopt32(argv, "");
 187         //argv += optind;
 188         argv++;
 189
 190         if (!*argv) {
 191                 /* Read from stdin, several numbers per line are accepted */
 192                 for (;;) {
 193                         char *numstr, *line;
 194                         line = xmalloc_fgetline(stdin);
 195                         if (!line)
 196                                 return EXIT_SUCCESS;
 197                         numstr = line;
 198                         for (;;) {
 199                                 char *end;
 200                                 numstr = skip_whitespace(numstr);
 201                                 if (!numstr[0])
 202                                         break;
 203                                 end = skip_non_whitespace(numstr);
 204                                 if (*end != '\0')
 205                                         *end++ = '\0';
 206                                 factorize_numstr(numstr);
 207                                 numstr = end;
 208                         }
 209                         free(line);
 210                 }
 211         }
 212
 213         do {
 214                 /* Leading spaces are ok (coreutils compat) */
 215                 factorize_numstr(skip_whitespace(*argv));
 216         } while (*++argv);
 217
 218         return EXIT_SUCCESS;
 219 }