src/base/low/cl_low_isqrt.cc

   1 // isqrt().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_low.h"
   8
   9
  10 // Implementation.
  11
  12 namespace cln {
  13
  14 // Zieht die Ganzzahl-Wurzel aus einer 32-Bit-Zahl und
  15 // liefert eine 16-Bit-Wurzel.
  16 // isqrt(x)
  17 // > uintL x : Radikand, >=0, <2^32
  18 // < uintL ergebnis : Wurzel, >=0, <2^16
  19 uintL isqrt (uintL x)
  20 {
  21   // Methode:
  22   // x=0 -> y=0, fertig.
  23   // y := 2^k als Anfangswert, wobei k>0, k<=16 mit 2^(2k-2) <= x < 2^(2k) sei.
  24   // y := floor((y + floor(x/y))/2) als nächster Wert,
  25   // solange z := floor(x/y) < y, setze y := floor((y+z)/2).
  26   // y ist fertig.
  27   // (Beweis:
  28   //  1. Die Folge der y ist streng monoton fallend.
  29   //  2. Stets gilt y >= floor(sqrt(x)) (denn für alle y>0 ist
  30   //     y + x/y >= 2*sqrt(x) und daher  floor((y + floor(x/y))/2) =
  31   //     floor(y/2 + x/(2*y)) >= floor(sqrt(x)) ).
  32   //  3. Am Schluß gilt x >= y^2.
  33   // )
  34      if (x==0) return 0; // x=0 -> y=0
  35      { var uintC k2; integerlength32(x,k2=); // 2^(k2-1) <= x < 2^k2
  36       {var uintC k1 = floor(k2-1,2); // k1 = k-1, k wie oben
  37        if (k1 < 16-1)
  38          // k < 16
  39          { var uintL y = (x >> (k1+2)) | bit(k1); // stets 2^(k-1) <= y < 2^k
  40            loop
  41              { var uintL z;
  42                divu_3216_1616(x,y, z=,); // Dividiere x/y (geht, da x/y < 2^(2k)/2^(k-1) = 2^(k+1) <= 2^16)
  43                if (z >= y) break;
  44                y = floor(z+y,2); // geht, da z+y < 2*y < 2^(k+1) <= 2^16
  45              }
  46            return y;
  47          }
  48          else
  49          // k = 16, Vorsicht!
  50          { var uintL x1 = high16(x);
  51            var uintL y = (x >> (16+1)) | bit(16-1); // stets 2^(k-1) <= y < 2^k
  52            loop
  53              { var uintL z;
  54                if (x1 >= y) break; // Division x/y ergäbe Überlauf -> z > y
  55                divu_3216_1616(x,y, z=,); // Dividiere x/y
  56                if (z >= y) break;
  57                y = floor(z+y,2);
  58              }
  59            return y;
  60          }
  61      }}
  62 }
  63
  64 #ifdef HAVE_LONGLONG
  65 uintL isqrt (uintQ x)
  66 {
  67   // As isqrt (uintL) above, but with 64-bit numbers.
  68      if (x==0) return 0; // x=0 -> y=0
  69      { var uintC k2; integerlength64(x,k2=); // 2^(k2-1) <= x < 2^k2
  70       {var uintC k1 = floor(k2-1,2); // k1 = k-1, k as above
  71        if (k1 < 32-1)
  72          // k < 32
  73          { var uintL y = (x >> (k1+2)) | bit(k1); // always 2^(k-1) <= y < 2^k
  74            loop
  75              { var uintL z;
  76                divu_6432_3232(high32(x),low32(x),y, z=,); // z := x/y (works, since x/y < 2^(2k)/2^(k-1) = 2^(k+1) <= 2^32)
  77                if (z >= y) break;
  78                y = floor(z+y,2); // geht, da z+y < 2*y < 2^(k+1) <= 2^32
  79              }
  80            return y;
  81          }
  82          else
  83          // k = 32, careful!
  84          { var uintL x1 = high32(x);
  85            var uintL y = (x >> (32+1)) | bit(32-1); // stets 2^(k-1) <= y < 2^k
  86            loop
  87              { var uintL z;
  88                if (x1 >= y) break; // division x/y would overflow -> z > y
  89                divu_6432_3232(high32(x),low32(x),y, z=,); // divide x/y
  90                if (z >= y) break;
  91                y = floor(z+y,2);
  92              }
  93            return y;
  94          }
  95      }}
  96 }
  97 #endif
  98
  99 }  // namespace cln