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"
   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 #define WIDE_BITS (unsigned)(sizeof(wide_t)*8)
  31 #define TOPMOST_WIDE_BIT ((wide_t)1 << (WIDE_BITS-1))
  32
  33 #if ULLONG_MAX == (UINT_MAX * UINT_MAX + 2 * UINT_MAX)
  34 /* "unsigned" is half as wide as ullong */
  35 typedef unsigned half_t;
  36 #define HALF_MAX UINT_MAX
  37 #define HALF_FMT ""
  38 #elif ULLONG_MAX == (ULONG_MAX * ULONG_MAX + 2 * ULONG_MAX)
  39 /* long is half as wide as ullong */
  40 typedef unsigned long half_t;
  41 #define HALF_MAX ULONG_MAX
  42 #define HALF_FMT "l"
  43 #else
  44 #error Cant find an integer type which is half as wide as ullong
  45 #endif
  46
  47 /* Returns such x that x+1 > sqrt(N) */
  48 static inline half_t isqrt(wide_t N)
  49 {
  50         wide_t mask_2bits;
  51         half_t x;
  52
  53 // Never called with N < 1
  54 //      if (N == 0)
  55 //              return 0;
  56
  57         /* First approximation of x+1 > sqrt(N) - all-ones, half as many bits:
  58          * 1xxxxx -> 111 (six bits to three)
  59          * 01xxxx -> 111
  60          * 001xxx -> 011
  61          * 0001xx -> 011 and so on.
  62          */
  63         x = HALF_MAX;
  64         mask_2bits = TOPMOST_WIDE_BIT | (TOPMOST_WIDE_BIT >> 1);
  65         while (!(N & mask_2bits)) {
  66                 x >>= 1;
  67                 mask_2bits >>= 2;
  68         }
  69         dbg("x:%"HALF_FMT"x", x);
  70
  71         for (;;) {
  72                 half_t y = (x + N/x) / 2;
  73                 dbg("y:%x y^2:%llx", y, (wide_t)y * y);
  74                 /*
  75                  * "real" y may be one bit wider: 0x100000000 and get truncated to 0.
  76                  * In this case, "real" y is > x. The first check below is for this case:
  77                  */
  78                 if (y == 0 || y >= x) {
  79                         dbg("isqrt(%llx)=%"HALF_FMT"x", N, x);
  80                         return x;
  81                 }
  82                 x = y;
  83         }
  84 }
  85
  86 static NOINLINE half_t isqrt_odd(wide_t N)
  87 {
  88         half_t s = isqrt(N);
  89         if (s && !(s & 1)) /* even? */
  90                 s--;
  91         return s;
  92 }
  93
  94 static NOINLINE void factorize(wide_t N)
  95 {
  96         half_t factor;
  97         half_t max_factor;
  98         unsigned count3;
  99         unsigned count5;
 100         unsigned count7;
 101
 102         if (N < 4)
 103                 goto end;
 104
 105         while (!(N & 1)) {
 106                 printf(" 2");
 107                 N >>= 1;
 108         }
 109
 110         max_factor = isqrt_odd(N);
 111         count3 = 3;
 112         count5 = 6;
 113         count7 = 9;
 114         factor = 3;
 115         for (;;) {
 116                 /* The division is the most costly part of the loop.
 117                  * On 64bit CPUs, takes at best 12 cycles, often ~20.
 118                  */
 119                 while ((N % factor) == 0) { /* not likely */
 120                         N = N / factor;
 121                         printf(" %"HALF_FMT"u", factor);
 122                         max_factor = isqrt_odd(N);
 123                 }
 124  next_factor:
 125                 if (factor >= max_factor)
 126                         break;
 127                 factor += 2;
 128                 /* Rudimentary wheel sieving: skip multiples of 3, 5 and 7:
 129                  * Every third odd number is divisible by three and thus isn't a prime:
 130                  * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47...
 131                  * ^ ^   ^  ^     ^  ^     ^  _     ^  ^     _  ^     ^  ^     ^
 132                  * (^ = primes, _ = would-be-primes-if-not-divisible-by-5)
 133                  */
 134                 count7--;
 135                 count5--;
 136                 count3--;
 137                 if (count3 && count5 && count7)
 138                         continue;
 139                 if (count3 == 0)
 140                         count3 = 3;
 141                 if (count5 == 0)
 142                         count5 = 5;
 143                 if (count7 == 0)
 144                         count7 = 7;
 145                 goto next_factor;
 146         }
 147  end:
 148         if (N > 1)
 149                 printf(" %llu", N);
 150         bb_putchar('\n');
 151 }
 152
 153 int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
 154 int factor_main(int argc UNUSED_PARAM, char **argv)
 155 {
 156         //// coreutils has undocumented option ---debug (three dashes)
 157         //getopt32(argv, "");
 158         //argv += optind;
 159         argv++;
 160
 161         if (!*argv)
 162                 //TODO: read from stdin
 163                 bb_show_usage();
 164
 165         do {
 166                 wide_t N;
 167                 const char *numstr;
 168
 169                 /* Coreutils compat */
 170                 numstr = skip_whitespace(*argv);
 171                 if (*numstr == '+')
 172                         numstr++;
 173
 174                 N = bb_strtoull(numstr, NULL, 10);
 175                 if (errno)
 176                         bb_show_usage();
 177                 printf("%llu:", N);
 178                 factorize(N);
 179         } while (*++argv);
 180
 181         return EXIT_SUCCESS;
 182 }