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 (3.2 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 #include "common_bufsiz.h"
  23
  24 #if 0
  25 # define dbg(...) bb_error_msg(__VA_ARGS__)
  26 #else
  27 # define dbg(...) ((void)0)
  28 #endif
  29
  30 typedef unsigned long long wide_t;
  31
  32 #if ULLONG_MAX == (UINT_MAX * UINT_MAX + 2 * UINT_MAX)
  33 /* "unsigned" is half as wide as ullong */
  34 typedef unsigned half_t;
  35 #define HALF_MAX UINT_MAX
  36 #define HALF_FMT ""
  37 #elif ULLONG_MAX == (ULONG_MAX * ULONG_MAX + 2 * ULONG_MAX)
  38 /* long is half as wide as ullong */
  39 typedef unsigned long half_t;
  40 #define HALF_MAX ULONG_MAX
  41 #define HALF_FMT "l"
  42 #else
  43 #error Cant find an integer type which is half as wide as ullong
  44 #endif
  45
  46 /* The trial divisor increment wheel.  Use it to skip over divisors that
  47  * are composites of 2, 3, 5, 7, or 11.
  48  * Larger wheels improve sieving only slightly, but quickly grow in size
  49  * (adding just one prime, 13, results in 5766 element sieve).
  50  */
  51 #define R(a,b,c,d,e,f,g,h,i,j,A,B,C,D,E,F,G,H,I,J,x) \
  52         (((uint64_t)(a<<0) | (b<<3) | (c<<6) | (d<<9) | (e<<12) | (f<<15) | (g<<18) | (h<<21) | (i<<24) | (j<<27)) << 1)  | \
  53         (((uint64_t)(A<<0) | (B<<3) | (C<<6) | (D<<9) | (E<<12) | (F<<15) | (G<<18) | (H<<21) | (I<<24) | (J<<27)) << 31) | \
  54         ((uint64_t)x << 61)
  55 #define P(a,b,c,d,e,f,g,h,i,j,A,B,C,D,E,F,G,H,I,J,x) \
  56         R(      (a/2),(b/2),(c/2),(d/2),(e/2),(f/2),(g/2),(h/2),(i/2),(j/2), \
  57                 (A/2),(B/2),(C/2),(D/2),(E/2),(F/2),(G/2),(H/2),(I/2),(J/2), \
  58                 (x/2) \
  59         )
  60 static const uint64_t packed_wheel[] = {
  61         /* 1, 2, */
  62         P( 2, 4, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6),
  63         P( 8, 4, 2, 4, 2, 4,14, 4, 6, 2,10, 2, 6, 6, 4, 2, 4, 6, 2,10, 2),
  64         P( 4, 2,12,10, 2, 4, 2, 4, 6, 2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4),
  65         P( 6, 8, 4, 2, 4, 6, 8, 6,10, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6),
  66         P( 4, 2, 6,10, 2,10, 2, 4, 2, 4, 6, 8, 4, 2, 4,12, 2, 6, 4, 2, 6),
  67         P( 4, 6,12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6,10, 2, 4, 6, 2),
  68         P( 6, 4, 2, 4, 2,10, 2,10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 4),
  69         P( 2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2,10),
  70         P( 2, 6, 4, 2, 4, 2,10, 2,10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2),
  71         P( 6, 4, 8, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4),
  72         P( 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2,10, 2, 6, 4, 6, 2, 6, 4, 2, 4),
  73         P( 6, 6, 8, 4, 2, 6,10, 8, 4, 2, 4, 2, 4, 8,10, 6, 2, 4, 8, 6, 6),
  74         P( 4, 2, 4, 6, 2, 6, 4, 6, 2,10, 2,10, 2, 4, 2, 4, 6, 2, 6, 4, 2),
  75         P( 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 8, 4, 6),
  76         P( 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2,10, 2,10, 2, 4, 2, 4, 6, 2),
  77         P(10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2),
  78         P( 4, 6, 2, 6, 6, 4, 6, 6, 2, 6, 6, 4, 2,10, 2,10, 2, 4, 2, 4, 6),
  79         P( 2, 6, 4, 2,10, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2,12, 6, 4),
  80         P( 6, 2, 4, 6, 2,12, 4, 2, 4, 8, 6, 4, 2, 4, 2,10, 2,10, 6, 2, 4),
  81         P( 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2,10, 6, 8, 6, 4, 2, 4, 8, 6),
  82         P( 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 2, 4, 2,10,12, 2, 4),
  83         P( 2,10, 2, 6, 4, 2, 4, 6, 6, 2,10, 2, 6, 4,14, 4, 2, 4, 2, 4, 8),
  84         P( 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4,12, 2,12),
  85 };
  86 #undef P
  87 #undef R
  88 #define WHEEL_START 5
  89 #define WHEEL_SIZE (5 + 24 * 20)
  90 #define square_count (((uint8_t*)&bb_common_bufsiz1)[0])
  91 #define wheel_tab    (((uint8_t*)&bb_common_bufsiz1) + 1)
  92 /*
  93  * Why, you ask?
  94  * plain byte array:
  95  *      function                old     new   delta
  96  *      wheel_tab                 -     485    +485
  97  * 3-bit-packed insanity:
  98  *      packed_wheel              -     184    +184
  99  *      factor_main             108     163     +55
 100  */
 101 static void unpack_wheel(void)
 102 {
 103         int i;
 104         uint8_t *p;
 105
 106         setup_common_bufsiz();
 107         wheel_tab[0] = 1;
 108         wheel_tab[1] = 2;
 109         p = &wheel_tab[2];
 110         for (i = 0; i < ARRAY_SIZE(packed_wheel); i++) {
 111                 uint64_t v = packed_wheel[i];
 112                 while ((v & 0xe) != 0) {
 113                         *p = v & 0xe;
 114                         //printf("%2u,", *p);
 115                         p++;
 116                         v >>= 3;
 117                 }
 118                 //printf("\n");
 119         }
 120 }
 121
 122 /* Prevent inlining, factorize() needs all help it can get with reducing register pressure */
 123 static NOINLINE void print_w(wide_t n)
 124 {
 125         unsigned rep = square_count;
 126         do
 127                 printf(" %llu", n);
 128         while (--rep != 0);
 129 }
 130 static NOINLINE void print_h(half_t n)
 131 {
 132         print_w(n);
 133 }
 134
 135 static void factorize(wide_t N);
 136
 137 static half_t isqrt_odd(wide_t N)
 138 {
 139         half_t s = isqrt(N);
 140         /* s^2 is <= N, (s+1)^2 > N */
 141
 142         /* If s^2 in fact is EQUAL to N, it's very lucky.
 143          * Examples:
 144          * factor 18446743988964486098 = 2 * 3037000493 * 3037000493
 145          * factor 18446743902517389507 = 3 * 2479700513 * 2479700513
 146          */
 147         if ((wide_t)s * s == N) {
 148                 /* factorize sqrt(N), printing each factor twice */
 149                 square_count *= 2;
 150                 factorize(s);
 151                 /* Let caller know we recursed */
 152                 return 0;
 153         }
 154
 155         /* Subtract 1 from even s, odd s won't change: */
 156         /* (doesnt work for zero, but we know that s != 0 here) */
 157         s = (s - 1) | 1;
 158         return s;
 159 }
 160
 161 static NOINLINE void factorize(wide_t N)
 162 {
 163         unsigned w;
 164         half_t factor;
 165         half_t max_factor;
 166
 167         if (N < 4)
 168                 goto end;
 169
 170         /* The code needs to be optimized for the case where
 171          * there are large prime factors. For example,
 172          * this is not hard:
 173          * 8262075252869367027 = 3 7 17 23 47 101 113 127 131 137 823
 174          *  (the largest divisor to test for largest factor 823
 175          *  is only ~sqrt(823) = 28, the entire factorization needs
 176          *  only ~33 trial divisions)
 177          * but this is:
 178          * 18446744073709551601 = 53 348051774975651917
 179          * the last factor requires testing up to
 180          * 589959129 - about 100 million iterations.
 181          * The slowest case (largest prime) for N < 2^64 is
 182          * factor 18446744073709551557 (0xffffffffffffffc5).
 183          */
 184         max_factor = isqrt_odd(N);
 185         if (!max_factor)
 186                 return; /* square was detected and recursively factored */
 187         factor = 2;
 188         w = 0;
 189         for (;;) {
 190                 half_t fw;
 191
 192                 /* The division is the most costly part of the loop.
 193                  * On 64bit CPUs, takes at best 12 cycles, often ~20.
 194                  */
 195                 while ((N % factor) == 0) { /* not likely */
 196                         N = N / factor;
 197                         print_h(factor);
 198                         max_factor = isqrt_odd(N);
 199                         if (!max_factor)
 200                                 return; /* square was detected */
 201                 }
 202                 if (factor >= max_factor)
 203                         break;
 204                 fw = factor + wheel_tab[w];
 205                 if (fw < factor)
 206                         break; /* overflow */
 207                 factor = fw;
 208                 w++;
 209                 if (w < WHEEL_SIZE)
 210                         continue;
 211                 w = WHEEL_START;
 212         }
 213  end:
 214         if (N > 1)
 215                 print_w(N);
 216         bb_putchar('\n');
 217 }
 218
 219 static void factorize_numstr(const char *numstr)
 220 {
 221         wide_t N;
 222
 223         /* Leading + is ok (coreutils compat) */
 224         if (*numstr == '+')
 225                 numstr++;
 226         N = bb_strtoull(numstr, NULL, 10);
 227         if (errno)
 228                 bb_show_usage();
 229         printf("%llu:", N);
 230         square_count = 1;
 231         factorize(N);
 232 }
 233
 234 int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
 235 int factor_main(int argc UNUSED_PARAM, char **argv)
 236 {
 237         unpack_wheel();
 238
 239         //// coreutils has undocumented option ---debug (three dashes)
 240         //getopt32(argv, "");
 241         //argv += optind;
 242         argv++;
 243
 244         if (!*argv) {
 245                 /* Read from stdin, several numbers per line are accepted */
 246                 for (;;) {
 247                         char *numstr, *line;
 248                         line = xmalloc_fgetline(stdin);
 249                         if (!line)
 250                                 return EXIT_SUCCESS;
 251                         numstr = line;
 252                         for (;;) {
 253                                 char *end;
 254                                 numstr = skip_whitespace(numstr);
 255                                 if (!numstr[0])
 256                                         break;
 257                                 end = skip_non_whitespace(numstr);
 258                                 if (*end != '\0')
 259                                         *end++ = '\0';
 260                                 factorize_numstr(numstr);
 261                                 numstr = end;
 262                         }
 263                         free(line);
 264                 }
 265         }
 266
 267         do {
 268                 /* Leading spaces are ok (coreutils compat) */
 269                 factorize_numstr(skip_whitespace(*argv));
 270         } while (*++argv);
 271
 272         return EXIT_SUCCESS;
 273 }