libbb/isqrt.c

   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) */
  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
  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