2 * Copyright (C) 2017 Denys Vlasenko <vda.linux@googlemail.com>
4 * Licensed under GPLv2, see file LICENSE in this source tree.
7 //config: bool "factor (2.6 kb)"
10 //config: factor factorizes integers
12 //applet:IF_FACTOR(APPLET(factor, BB_DIR_USR_BIN, BB_SUID_DROP))
14 //kbuild:lib-$(CONFIG_FACTOR) += factor.o
16 //usage:#define factor_trivial_usage
17 //usage: "[NUMBER]..."
18 //usage:#define factor_full_usage "\n\n"
19 //usage: "Print prime factors"
24 # define dbg(...) bb_error_msg(__VA_ARGS__)
26 # define dbg(...) ((void)0)
29 typedef unsigned long long wide_t;
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
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
42 #error Cant find an integer type which is half as wide as ullong
45 static half_t isqrt_odd(wide_t N)
48 /* Subtract 1 from even s, odd s won't change: */
49 /* (doesnt work for zero, but we know that s != 0 here) */
54 static NOINLINE void factorize(wide_t N)
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
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,
88 /* The code needs to be optimized for the case where
89 * there are large prime factors. For example,
91 * 8262075252869367027 = 3 7 17 23 47 101 113 127 131 137 823
92 * (the largest factor to test is only ~sqrt(823) = 28)
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).
100 max_factor = isqrt_odd(N);
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)
115 /* The division is the most costly part of the loop.
116 * On 64bit CPUs, takes at best 12 cycles, often ~20.
118 while ((N % factor) == 0) { /* not likely */
120 printf(" %"HALF_FMT"u", factor);
121 max_factor = isqrt_odd(N);
124 if (factor >= max_factor)
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.
137 // if (count3 && count5 && count7)
139 sieve_word += INCREMENT_EACH;
140 if (!(sieve_word & MULTIPLE_DETECTED))
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.
150 if (sieve_word & MULTIPLE_OF_3)
151 sieve_word -= SHIFT_3 * 3;
154 if (sieve_word & MULTIPLE_OF_5)
155 sieve_word -= SHIFT_5 * 5;
158 if (sieve_word & MULTIPLE_OF_7)
159 sieve_word -= SHIFT_7 * 7;
168 static void factorize_numstr(const char *numstr)
172 /* Leading + is ok (coreutils compat) */
175 N = bb_strtoull(numstr, NULL, 10);
182 int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
183 int factor_main(int argc UNUSED_PARAM, char **argv)
185 //// coreutils has undocumented option ---debug (three dashes)
186 //getopt32(argv, "");
191 /* Read from stdin, several numbers per line are accepted */
194 line = xmalloc_fgetline(stdin);
200 numstr = skip_whitespace(numstr);
203 end = skip_non_whitespace(numstr);
206 factorize_numstr(numstr);
214 /* Leading spaces are ok (coreutils compat) */
215 factorize_numstr(skip_whitespace(*argv));