* 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"
#endif
typedef unsigned long long wide_t;
-#define WIDE_BITS (unsigned)(sizeof(wide_t)*8)
-#define TOPMOST_WIDE_BIT ((wide_t)1 << (WIDE_BITS-1))
#if ULLONG_MAX == (UINT_MAX * UINT_MAX + 2 * UINT_MAX)
/* "unsigned" is half as wide as ullong */
#error Cant find an integer type which is half as wide as ullong
#endif
-/* Returns such x that x+1 > sqrt(N) */
-static inline half_t isqrt(wide_t N)
-{
- half_t x;
-
- // Never called with N < 1
- //if (N == 0)
- // return 0;
-
- x = HALF_MAX;
- /* First approximation of x+1 > sqrt(N) - all-ones, half as many bits:
- * 1xxxxx -> 111 (six bits to three)
- * 01xxxx -> 111
- * 001xxx -> 011
- * 0001xx -> 011 and so on.
- */
- // It is actually not performance-critical at all.
- // Can simply start Newton loop with very conservative x=0xffffffff.
- //wide_t mask_2bits;
- //mask_2bits = TOPMOST_WIDE_BIT | (TOPMOST_WIDE_BIT >> 1);
- //while (!(N & mask_2bits)) {
- // x >>= 1;
- // mask_2bits >>= 2;
- //}
- dbg("x:%"HALF_FMT"x", x);
-
- for (;;) {
- half_t y = (x + N/x) / 2;
- dbg("y:%x y^2:%llx", y, (wide_t)y * y);
- /*
- * "real" y may be one bit wider: 0x100000000 and get truncated to 0.
- * In this case, "real" y is > x. The first check below is for this case:
- */
- if (y == 0 || y >= x) {
- dbg("isqrt(%llx)=%"HALF_FMT"x", N, x);
- return x;
- }
- x = y;
- }
-}
-
-static NOINLINE half_t isqrt_odd(wide_t N)
+static half_t isqrt_odd(wide_t N)
{
half_t s = isqrt(N);
/* Subtract 1 from even s, odd s won't change: */
// unsigned count3;
// unsigned count5;
// unsigned count7;
- // ^^^^^^^^^^^^^^^ commented-out simple siving code (easier to grasp).
+ // ^^^^^^^^^^^^^^^ 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
MULTIPLE_OF_3 = 1 << 2,
MULTIPLE_OF_5 = 1 << 6,
MULTIPLE_OF_7 = 1 << 11,
- MULTIPLE_3_5_7 = MULTIPLE_OF_3 | MULTIPLE_OF_5 | MULTIPLE_OF_7,
+ MULTIPLE_DETECTED = MULTIPLE_OF_3 | MULTIPLE_OF_5 | MULTIPLE_OF_7,
};
unsigned sieve_word;
* 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;
for (;;) {
// if (count3 && count5 && count7)
// continue;
sieve_word += INCREMENT_EACH;
- if (!(sieve_word & MULTIPLE_3_5_7))
+ if (!(sieve_word & MULTIPLE_DETECTED))
continue;
/*
* "factor" is multiple of 3 33% of the time (count3 reached 0),
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)
{
-#if 0 /* test isqrt() correctness */
- wide_t n = argv[1] ? bb_strtoull(argv[1], NULL, 0) : time(NULL);
- for (;;) {
- half_t h;
- if (--n == 0)
- --n;
- h = isqrt(n);
- if (!(n & 0xff))
- printf("isqrt(%llx)=%"HALF_FMT"x\n", n, h);
- if ((wide_t)h * h > n)
- return 1;
- h++;
- if (h != 0 && (wide_t)h * h <= n)
- return 1;
- }
-#endif
-
//// coreutils has undocumented option ---debug (three dashes)
//getopt32(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 {
- wide_t N;
- const char *numstr;
-
- /* Coreutils compat */
- numstr = skip_whitespace(*argv);
- if (*numstr == '+')
- numstr++;
-
- N = bb_strtoull(numstr, NULL, 10);
- if (errno)
- bb_show_usage();
- printf("%llu:", N);
- factorize(N);
+ /* Leading spaces are ok (coreutils compat) */
+ factorize_numstr(skip_whitespace(*argv));
} while (*++argv);
return EXIT_SUCCESS;
projects / oweals / busybox.git / blobdiff