* Licensed under GPLv2, see file LICENSE in this source tree.
*/
//config:config FACTOR
-//config: bool "factor"
+//config: bool "factor (2.6 kb)"
//config: default y
//config: help
-//config: factor factorizes integers
+//config: factor factorizes integers
//applet:IF_FACTOR(APPLET(factor, BB_DIR_USR_BIN, BB_SUID_DROP))
//kbuild:lib-$(CONFIG_FACTOR) += factor.o
//usage:#define factor_trivial_usage
-//usage: "NUMBER..."
+//usage: "[NUMBER]..."
//usage:#define factor_full_usage "\n\n"
//usage: "Print prime factors"
#include "libbb.h"
+#if 0
+# define dbg(...) bb_error_msg(__VA_ARGS__)
+#else
+# define dbg(...) ((void)0)
+#endif
+
+typedef unsigned long long wide_t;
+
#if ULLONG_MAX == (UINT_MAX * UINT_MAX + 2 * UINT_MAX)
/* "unsigned" is half as wide as ullong */
typedef unsigned half_t;
+#define HALF_MAX UINT_MAX
#define HALF_FMT ""
#elif ULLONG_MAX == (ULONG_MAX * ULONG_MAX + 2 * ULONG_MAX)
/* long is half as wide as ullong */
typedef unsigned long half_t;
+#define HALF_MAX ULONG_MAX
#define HALF_FMT "l"
#else
#error Cant find an integer type which is half as wide as ullong
#endif
-static void factorize(const char *numstr)
+static half_t isqrt_odd(wide_t N)
{
- unsigned long long N, factor2;
- half_t factor;
- unsigned count3;
-
- /* Coreutils compat */
- numstr = skip_whitespace(numstr);
- if (*numstr == '+')
- numstr++;
-
- N = bb_strtoull(numstr, NULL, 10);
- if (errno)
- bb_show_usage();
+ half_t s = isqrt(N);
+ /* Subtract 1 from even s, odd s won't change: */
+ /* (doesnt work for zero, but we know that s != 0 here) */
+ s = (s - 1) | 1;
+ return s;
+}
- printf("%llu:", N);
+static NOINLINE void factorize(wide_t N)
+{
+ half_t factor;
+ half_t max_factor;
+ // unsigned count3;
+ // unsigned count5;
+ // unsigned count7;
+ // ^^^^^^^^^^^^^^^ commented-out simple sieving code (easier to grasp).
+ // Faster sieving, using one word for potentially up to 6 counters:
+ // count upwards in each mask, counter "triggers" when it sets its mask to "100[0]..."
+ // 10987654321098765432109876543210 - bits 31-0 in 32-bit word
+ // 17777713333311111777775555333 - bit masks for counters for primes 3,5,7,11,13,17
+ // 100000100001000010001001 - value for adding 1 to each mask
+ // 10000010000010000100001000100 - value for checking that any mask reached msb
+ enum {
+ SHIFT_3 = 1 << 0,
+ SHIFT_5 = 1 << 3,
+ SHIFT_7 = 1 << 7,
+ INCREMENT_EACH = SHIFT_3 | SHIFT_5 | SHIFT_7,
+ MULTIPLE_OF_3 = 1 << 2,
+ MULTIPLE_OF_5 = 1 << 6,
+ MULTIPLE_OF_7 = 1 << 11,
+ MULTIPLE_DETECTED = MULTIPLE_OF_3 | MULTIPLE_OF_5 | MULTIPLE_OF_7,
+ };
+ unsigned sieve_word;
if (N < 4)
goto end;
+
while (!(N & 1)) {
printf(" 2");
N >>= 1;
}
- count3 = 3;
+ /* The code needs to be optimized for the case where
+ * there are large prime factors. For example,
+ * this is not hard:
+ * 8262075252869367027 = 3 7 17 23 47 101 113 127 131 137 823
+ * (the largest factor to test is only ~sqrt(823) = 28)
+ * but this is:
+ * 18446744073709551601 = 53 348051774975651917
+ * the last factor requires testing up to
+ * 589959129 - about 100 million iterations.
+ * The slowest case (largest prime) for N < 2^64 is
+ * factor 18446744073709551557 (0xffffffffffffffc5).
+ */
+ max_factor = isqrt_odd(N);
+ // count3 = 3;
+ // count5 = 6;
+ // count7 = 9;
+ sieve_word = 0
+ /* initial count for SHIFT_n is (n-1)/2*3: */
+ + (MULTIPLE_OF_3 - 3 * SHIFT_3)
+ + (MULTIPLE_OF_5 - 6 * SHIFT_5)
+ + (MULTIPLE_OF_7 - 9 * SHIFT_7)
+ //+ (MULTIPLE_OF_11 - 15 * SHIFT_11)
+ //+ (MULTIPLE_OF_13 - 18 * SHIFT_13)
+ //+ (MULTIPLE_OF_17 - 24 * SHIFT_17)
+ ;
factor = 3;
- factor2 = 3 * 3;
for (;;) {
- unsigned long long nfactor2;
-
- while ((N % factor) == 0) {
+ /* The division is the most costly part of the loop.
+ * On 64bit CPUs, takes at best 12 cycles, often ~20.
+ */
+ while ((N % factor) == 0) { /* not likely */
N = N / factor;
- printf(" %u"HALF_FMT"", factor);
+ printf(" %"HALF_FMT"u", factor);
+ max_factor = isqrt_odd(N);
}
next_factor:
- /* (f + 2)^2 = f^2 + 4*f + 4 = f^2 + 4*(f+1) */
- nfactor2 = factor2 + 4 * (factor + 1);
- if (nfactor2 < factor2) /* overflow? */
- break;
- factor2 = nfactor2;
- if (factor2 > N)
+ if (factor >= max_factor)
break;
factor += 2;
- /* Rudimentary wheel sieving: skip multiples of 3:
+ /* Rudimentary wheel sieving: skip multiples of 3, 5 and 7:
* Every third odd number is divisible by three and thus isn't a prime:
- * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37...
- * ^ ^ ^ ^ ^ ^ ^ _ ^ ^ _ ^ (^primes)
+ * 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47...
+ * ^ ^ ^ ^ ^ ^ ^ _ ^ ^ _ ^ ^ ^ ^
+ * (^ = primes, _ = would-be-primes-if-not-divisible-by-5)
+ * The numbers with space under them are excluded by sieve 3.
*/
- --count3;
- if (count3 == 0) {
- count3 = 3;
- goto next_factor;
- }
+ // count7--;
+ // count5--;
+ // count3--;
+ // if (count3 && count5 && count7)
+ // continue;
+ sieve_word += INCREMENT_EACH;
+ if (!(sieve_word & MULTIPLE_DETECTED))
+ continue;
+ /*
+ * "factor" is multiple of 3 33% of the time (count3 reached 0),
+ * else, multiple of 5 13% of the time,
+ * else, multiple of 7 7.6% of the time.
+ * Cumulatively, with 3,5,7 sieving we are here 54.3% of the time.
+ */
+ // if (count3 == 0)
+ // count3 = 3;
+ if (sieve_word & MULTIPLE_OF_3)
+ sieve_word -= SHIFT_3 * 3;
+ // if (count5 == 0)
+ // count5 = 5;
+ if (sieve_word & MULTIPLE_OF_5)
+ sieve_word -= SHIFT_5 * 5;
+ // if (count7 == 0)
+ // count7 = 7;
+ if (sieve_word & MULTIPLE_OF_7)
+ sieve_word -= SHIFT_7 * 7;
+ goto next_factor;
}
end:
if (N > 1)
bb_putchar('\n');
}
+static void factorize_numstr(const char *numstr)
+{
+ wide_t N;
+
+ /* Leading + is ok (coreutils compat) */
+ if (*numstr == '+')
+ numstr++;
+ N = bb_strtoull(numstr, NULL, 10);
+ if (errno)
+ bb_show_usage();
+ printf("%llu:", N);
+ factorize(N);
+}
+
int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
int factor_main(int argc UNUSED_PARAM, char **argv)
{
//argv += optind;
argv++;
- if (!*argv)
- //TODO: read from stdin
- bb_show_usage();
+ if (!*argv) {
+ /* Read from stdin, several numbers per line are accepted */
+ for (;;) {
+ char *numstr, *line;
+ line = xmalloc_fgetline(stdin);
+ if (!line)
+ return EXIT_SUCCESS;
+ numstr = line;
+ for (;;) {
+ char *end;
+ numstr = skip_whitespace(numstr);
+ if (!numstr[0])
+ break;
+ end = skip_non_whitespace(numstr);
+ if (*end != '\0')
+ *end++ = '\0';
+ factorize_numstr(numstr);
+ numstr = end;
+ }
+ free(line);
+ }
+ }
do {
- factorize(*argv);
+ /* Leading spaces are ok (coreutils compat) */
+ factorize_numstr(skip_whitespace(*argv));
} while (*++argv);
- fflush_stdout_and_exit(EXIT_SUCCESS);
+ return EXIT_SUCCESS;
}