3 files changed, 62 insertions, 5 deletions

diff --git a/coreutils/factor.c b/coreutils/factor.c
index 90e9ed761..4df523f28 100644
--- a/coreutils/factor.c
+++ b/coreutils/factor.c
@@ -136,7 +136,7 @@ static void factorize(wide_t N);
 
 static half_t isqrt_odd(wide_t N)
 {
-        half_t s = isqrt(N);
+        half_t s = isqrt_ull(N);
         /* s^2 is <= N, (s+1)^2 > N */
 
         /* If s^2 in fact is EQUAL to N, it's very lucky.
diff --git a/include/libbb.h b/include/libbb.h
index aacff42ee..7cca925b9 100644
--- a/include/libbb.h
+++ b/include/libbb.h
@@ -408,7 +408,12 @@ unsigned FAST_FUNC bb_popcnt_long(unsigned long m);
 #define bb_popcnt_long(m) bb_popcnt_32(m)
 #endif
 
-unsigned long FAST_FUNC isqrt(unsigned long long N);
+unsigned FAST_FUNC isqrt(unsigned long N);
+#if LLONG_MAX > LONG_MAX
+unsigned long FAST_FUNC isqrt_ull(unsigned long long N);
+#else
+# define isqrt_ull(N) isqrt(N)
+#endif
 
 unsigned long long monotonic_ns(void) FAST_FUNC;
 unsigned long long monotonic_us(void) FAST_FUNC;
diff --git a/libbb/isqrt.c b/libbb/isqrt.c
index 817b7d036..507c7d637 100644
--- a/libbb/isqrt.c
+++ b/libbb/isqrt.c
@@ -9,15 +9,42 @@
 #ifndef ISQRT_TEST
 # include "libbb.h"
 #else
-/* gcc -DISQRT_TEST -Wall -O2 isqrt.c -oisqrt && ./isqrt $((RANDOM*RANDOM)) */
+// gcc -DISQRT_TEST -Wall -Os -ffunction-sections isqrt.c -S -fverbose-asm
+// gcc -DISQRT_TEST -Wall -Os -ffunction-sections isqrt.c -c
+// gcc -DISQRT_TEST -Wall -Os -ffunction-sections isqrt.c -oisqrt && ./isqrt $((RANDOM*RANDOM))
 # include <stdlib.h>
 # include <stdio.h>
 # include <time.h>
+# include <limits.h>
 # define FAST_FUNC /* nothing */
 #endif
 
 /* Returns such x that x+1 > sqrt(N) */
-unsigned long FAST_FUNC isqrt(unsigned long long N)
+
+/* Stolen from kernel code */
+unsigned FAST_FUNC isqrt(unsigned long x)
+{
+        unsigned long y = 0;
+        unsigned long m = 1UL << (8*sizeof(x) - 2);
+        // ^^^^^ should be m = 1UL << ((fls(x) - 1) & ~1UL);
+        // but calculating fls() would take time and code
+
+        do {
+                unsigned long b = y + m;
+                y >>= 1;
+                if (x >= b) {
+                        x -= b;
+                        y += m;
+                }
+                m >>= 2;
+        } while (m);
+        return y;
+}
+
+#if LLONG_MAX > LONG_MAX
+# if defined(i386)
+/* Smaller code on this register-starved arch */
+unsigned long FAST_FUNC isqrt_ull(unsigned long long N)
 {
         unsigned long x;
         unsigned shift;
@@ -33,15 +60,40 @@ unsigned long FAST_FUNC isqrt(unsigned long long N)
         } while ((int)shift >= 0);
         return x;
 }
+# else
+/* Stolen from kernel code */
+unsigned long FAST_FUNC isqrt_ull(unsigned long long x)
+{
+        unsigned long long y = 0;
+        unsigned long long m = 1ULL << (8*sizeof(x) - 2);
+        // ^^^^^ should be m = 1ULL << ((fls(x) - 1) & ~1ULL);
+        // but calculating fls() would take time and code
+
+        do {
+                unsigned long long b = y + m;
+                y >>= 1;
+                if (x >= b) {
+                        x -= b;
+                        y += m;
+                }
+                m >>= 2;
+        } while (m);
+        return y;
+}
+# endif
+#endif /* wide version needed */
 
 #ifdef ISQRT_TEST
+# if LLONG_MAX == LONG_MAX
+#  define isqrt_ull(N) isqrt(N)
+# endif
 int main(int argc, char **argv)
 {
         unsigned long long n = argv[1] ? strtoull(argv[1], NULL, 0) : time(NULL);
         for (;;) {
                 unsigned long h;
                 n--;
-                h = isqrt(n);
+                h = isqrt_ull(n);
                 if (!(n & 0xffff))
                         printf("isqrt(%llx)=%lx\n", n, h);
                 if ((unsigned long long)h * h > n) {