1 /*
2 * Copyright (C) 2017 Denys Vlasenko <vda.linux@googlemail.com>
3 *
4 * Licensed under GPLv2, see file LICENSE in this source tree.
5 */
6
7 //kbuild:lib-y += isqrt.o
8
9 #ifndef ISQRT_TEST
10 # include "libbb.h"
11 #else
12 /* gcc -DISQRT_TEST -Wall -O2 isqrt.c -oisqrt && ./isqrt $((RANDOM*RANDOM)) */
13 # include <stdlib.h>
14 # include <stdio.h>
15 # include <time.h>
16 # define FAST_FUNC /* nothing */
17 #endif
18
19 /* Returns such x that x+1 > sqrt(N) */
isqrt(unsigned long long N)20 unsigned long FAST_FUNC isqrt(unsigned long long N)
21 {
22 unsigned long x;
23 unsigned shift;
24 #define LL_WIDTH_BITS (unsigned)(sizeof(N)*8)
25
26 shift = LL_WIDTH_BITS - 2;
27 x = 0;
28 do {
29 x = (x << 1) + 1;
30 if ((unsigned long long)x * x > (N >> shift))
31 x--; /* whoops, that +1 was too much */
32 shift -= 2;
33 } while ((int)shift >= 0);
34 return x;
35 }
36
37 #ifdef ISQRT_TEST
main(int argc,char ** argv)38 int main(int argc, char **argv)
39 {
40 unsigned long long n = argv[1] ? strtoull(argv[1], NULL, 0) : time(NULL);
41 for (;;) {
42 unsigned long h;
43 n--;
44 h = isqrt(n);
45 if (!(n & 0xffff))
46 printf("isqrt(%llx)=%lx\n", n, h);
47 if ((unsigned long long)h * h > n) {
48 printf("BAD1: isqrt(%llx)=%lx\n", n, h);
49 return 1;
50 }
51 h++;
52 if ((unsigned long long)h * h != 0 /* this can overflow to 0 - not a bug */
53 && (unsigned long long)h * h <= n)
54 {
55 printf("BAD2: isqrt(%llx)=%lx\n", n, h);
56 return 1;
57 }
58 }
59 }
60 #endif
61