aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorDenys Vlasenko <vda.linux@googlemail.com>2017-04-10 00:24:16 +0200
committerDenys Vlasenko <vda.linux@googlemail.com>2017-04-10 00:24:16 +0200
commitc6476dca5484a9413d20dcbf887cc78055e6965f (patch)
tree96c1328084f6658e9d77b1abad056d0bb4feefc6
parentc804d4ec5cb222c842644bb99d9b077f5c6576f2 (diff)
downloadbusybox-c6476dca5484a9413d20dcbf887cc78055e6965f.tar.gz
factor: simpler isqrt
function old new delta isqrt_odd 102 87 -15 Signed-off-by: Denys Vlasenko <vda.linux@googlemail.com>
-rw-r--r--coreutils/factor.c63
1 files changed, 39 insertions, 24 deletions
diff --git a/coreutils/factor.c b/coreutils/factor.c
index 85284aa27..48833e103 100644
--- a/coreutils/factor.c
+++ b/coreutils/factor.c
@@ -47,33 +47,37 @@ typedef unsigned long half_t;
/* Returns such x that x+1 > sqrt(N) */
static inline half_t isqrt(wide_t N)
{
- wide_t x;
- unsigned c;
+ wide_t mask_2bits;
+ half_t x;
// Never called with N < 1
// if (N == 0)
// return 0;
-//
- /* Count leading zeros */
- c = 0;
- while (!(N & TOPMOST_WIDE_BIT)) {
- c++;
- N <<= 1;
+
+ /* 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;
}
- N >>= c;
+ dbg("x:%"HALF_FMT"x", x);
- /* Make x > sqrt(n) */
- x = (wide_t)1 << ((WIDE_BITS + 1 - c) / 2);
- dbg("x:%llx", x);
for (;;) {
- wide_t y = (x + N/x) / 2;
- dbg("y:%llx y^2:%llx (y+1)^2:%llx]", y, y*y, (y+1)*(y+1));
- if (y >= x) {
- /* Handle degenerate case N = 0xffffffffff...fffffff */
- if (y == (wide_t)HALF_MAX + 1)
- y--;
- dbg("isqrt(%llx)=%llx"HALF_FMT, N, y);
- return y;
+ 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;
}
@@ -92,6 +96,8 @@ static NOINLINE void factorize(wide_t N)
half_t factor;
half_t max_factor;
unsigned count3;
+ unsigned count5;
+ unsigned count7;
if (N < 4)
goto end;
@@ -103,6 +109,8 @@ static NOINLINE void factorize(wide_t N)
max_factor = isqrt_odd(N);
count3 = 3;
+ count5 = 6;
+ count7 = 9;
factor = 3;
for (;;) {
/* The division is the most costly part of the loop.
@@ -123,11 +131,18 @@ static NOINLINE void factorize(wide_t N)
* ^ ^ ^ ^ ^ ^ ^ _ ^ ^ _ ^ ^ ^ ^
* (^ = primes, _ = would-be-primes-if-not-divisible-by-5)
*/
- --count3;
- if (count3 == 0) {
+ count7--;
+ count5--;
+ count3--;
+ if (count3 && count5 && count7)
+ continue;
+ if (count3 == 0)
count3 = 3;
- goto next_factor;
- }
+ if (count5 == 0)
+ count5 = 5;
+ if (count7 == 0)
+ count7 = 7;
+ goto next_factor;
}
end:
if (N > 1)