3
0

factor.c 5.7 KB

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