X-Git-Url: https://git.librecmc.org/?a=blobdiff_plain;f=coreutils%2Ffactor.c;h=467e23a4d81aebbda0fcba7cab6f31d05f14c0d8;hb=feb79e8742eb3cef211804dadcc7f3ddfd154c72;hp=48833e10361a389cdf84632c427dd78a54635379;hpb=c6476dca5484a9413d20dcbf887cc78055e6965f;p=oweals%2Fbusybox.git diff --git a/coreutils/factor.c b/coreutils/factor.c index 48833e103..467e23a4d 100644 --- a/coreutils/factor.c +++ b/coreutils/factor.c @@ -4,17 +4,17 @@ * 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" @@ -27,8 +27,6 @@ #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 */ @@ -44,50 +42,12 @@ typedef unsigned long half_t; #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) -{ - wide_t mask_2bits; - half_t x; - -// Never called with N < 1 -// if (N == 0) -// return 0; - - /* First approximation x > sqrt(N) - half as many bits: - * 1xxxxx -> 111 (six bits to three) - * 01xxxx -> 111 - * 001xxx -> 011 - * 0001xx -> 011 and so on. - */ - x = HALF_MAX; - mask_2bits = TOPMOST_WIDE_BIT | (TOPMOST_WIDE_BIT >> 1); - while (mask_2bits && !(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); - if (s && !(s & 1)) /* even? */ - s--; + /* 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; } @@ -95,9 +55,27 @@ static NOINLINE void factorize(wide_t N) { half_t factor; half_t max_factor; - unsigned count3; - unsigned count5; - unsigned count7; + // 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; @@ -107,10 +85,31 @@ static NOINLINE void factorize(wide_t N) N >>= 1; } + /* 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; + // 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 (;;) { /* The division is the most costly part of the loop. @@ -118,30 +117,46 @@ static NOINLINE void factorize(wide_t N) */ 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: 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 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. */ - count7--; - count5--; - count3--; - if (count3 && count5 && count7) + // count7--; + // count5--; + // count3--; + // if (count3 && count5 && count7) + // continue; + sieve_word += INCREMENT_EACH; + if (!(sieve_word & MULTIPLE_DETECTED)) continue; - if (count3 == 0) - count3 = 3; - if (count5 == 0) - count5 = 5; - if (count7 == 0) - count7 = 7; + /* + * "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: @@ -150,6 +165,20 @@ static NOINLINE void factorize(wide_t N) 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) { @@ -158,24 +187,32 @@ 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 { - 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;