Method sqrt

ext/src/llvm/APInt.cpp:1097–1169

Square Root - this method computes and returns the square root of "this". Three mechanisms are used for computation. For small values (<= 5 bits), a table lookup is done. This gets some performance for common cases. For values using less than 52 bits, the value is converted to double and then the libc sqrt function is called. The result is rounded and then converted back to a uint64_t which is the

// back to a uint64_t which is then used to construct the result. Finally,
// the Babylonian method for computing square roots is used.
APInt APInt::sqrt() const {

  // Determine the magnitude of the value.
  unsigned magnitude = getActiveBits();

  // Use a fast table for some small values. This also gets rid of some
  // rounding errors in libc sqrt for small values.
  if (magnitude <= 5) {
    static const uint8_t results[32] = {
      /*     0 */ 0,
      /*  1- 2 */ 1, 1,
      /*  3- 6 */ 2, 2, 2, 2,
      /*  7-12 */ 3, 3, 3, 3, 3, 3,
      /* 13-20 */ 4, 4, 4, 4, 4, 4, 4, 4,
      /* 21-30 */ 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
      /*    31 */ 6
    };
    return APInt(BitWidth, results[ (isSingleWord() ? U.VAL : U.pVal[0]) ]);
  }

  // If the magnitude of the value fits in less than 52 bits (the precision of
  // an IEEE double precision floating point value), then we can use the
  // libc sqrt function which will probably use a hardware sqrt computation.
  // This should be faster than the algorithm below.
  if (magnitude < 52) {
    return APInt(BitWidth,
                 uint64_t(::round(::sqrt(double(isSingleWord() ? U.VAL
                                                               : U.pVal[0])))));
  }

  // Okay, all the short cuts are exhausted. We must compute it. The following
  // is a classical Babylonian method for computing the square root. This code
  // was adapted to APInt from a wikipedia article on such computations.
  // See http://www.wikipedia.org/ and go to the page named
  // Calculate_an_integer_square_root.
  unsigned nbits = BitWidth, i = 4;
  APInt testy(BitWidth, 16);
  APInt x_old(BitWidth, 1);
  APInt x_new(BitWidth, 0);
  APInt two(BitWidth, 2);

  // Select a good starting value using binary logarithms.
  for (;; i += 2, testy = testy.shl(2))
    if (i >= nbits || this->ule(testy)) {
      x_old = x_old.shl(i / 2);
      break;
    }

  // Use the Babylonian method to arrive at the integer square root:
  for (;;) {
    x_new = (this->udiv(x_old) + x_old).udiv(two);
    if (x_old.ule(x_new))
      break;
    x_old = x_new;
  }

  // Make sure we return the closest approximation
  // NOTE: The rounding calculation below is correct. It will produce an