calgorithmsquare-root

Square root source code found on the net


I'm looking for some algorithm for square root calculation and found this source file. I would like to try to replicate it because it seems quite simple but I can not relate it to some known algorithm (Newton, Babylon ...). Can you tell me the name?

int sqrt(int num) {
    int op = num;
    int res = 0;
    int one = 1 << 30; // The second-to-top bit is set: 1L<<30 for long

    // "one" starts at the highest power of four <= the argument.
    while (one > op)
        one >>= 2;

    while (one != 0) {
        if (op >= res + one) {
            op -= res + one;
            res += 2 * one;
        }
        res >>= 1;
        one >>= 2;
    }
    return res;
}

Solution

  • As @Eugene Sh. references, this is the classic "digit-by-digit" method done to calculate square root. Learned it in base 10 when such things were taught in primary school.

    OP's code fails select numbers too. sqrt(1073741824) --> -1 rather than expected 32768. 1073741824 == 0x40000000. Further, it fails most (all?) values this and greater. Of course OP's sqrt(some_negative) is a problem too.

    Candidate alternative: also here

    unsigned isqrt(unsigned num) {
      unsigned res = 0;
      
      // The second-to-top bit is set: 1 << 30 for 32 bits
      // Needs work to run on unusual platforms where `unsigned` has padding or odd bit width.
      unsigned bit = 1u << (sizeof(num) * CHAR_BIT - 2); 
    
      // "bit" starts at the highest power of four <= the argument.
      while (bit > num) {
        bit >>= 2;
      }
    
      while (bit > 0) {
        if (num >= res + bit) {
          num -= res + bit;
          res = (res >> 1) + bit;  // Key difference between this and OP's code
        } else {
          res >>= 1;
        }
        bit >>= 2;
      }
    
      return res;
    }
    

    Portability update. The greatest power of 4 is needed.

    #include <limits.h>
    // greatest power of 4 <= a power-of-2 minus 1
    #define POW4_LE_POW2M1(n) (  ((n)/2 + 1) >> ((n)%3==0)  )
    
    unsigned bit = POW4_LE_POW2M1(UINT_MAX);