```   1/*
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
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
29typedef 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 */
33typedef 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 */
38typedef 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
45static 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
54static 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
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
168static 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
182int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
183int 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}
220```