Help to understand algorithm for finding square root?

Question

Below there is an excerpt of an algorithm to compute square roots. I present it as I found.

Step 1: Group the number in "twos" from the decimal place. (If you have a number with an odd number of digits, the group to the far left will only be a group of 1 digit.)

Step 2: Start with the first group of two (the group on the left). This group may be a group of only one number if your number has an odd number of digits. Find the greatest square less than or equal to that group of digits and its square root will be your first approximation of the entire square root.

Step 3: Subtract the current approximation squared and bring down the next group of numbers behind it. This is your next number to work with.

Step 4: Double the current approximation of the root.

Step 5: Find the "ones" digit of the doubled number that would result in a number which divides into the number you are currently working with- with the smallest possible remainder. This is the next number in your approximation of the square root.

Step 6: Multiply the "ones" digit by the doubled number plus the "ones" digit. Subtract this number from the number you are currently working with and bring down the next group of numbers behind it. This is your next group of numbers to work with.

Step 7: Repeat steps 4 through 6 until you get an approximation with an acceptable number of significant digits.

I understand everything but Step 5. In the picture above, for instance, how do they get that $6$ in order to form $26$?

I had understood we should find the digit $\bullet$ in $2\bullet$ such that $154/2\bullet$ would produce the least remainder. But if that were the case, one should have $\bullet=9$ since $154/29$ has remainder $0$.

Furtheremore, what if there were more than one digit $\bullet$ which lead to the same remainder? Which would I pick?

Thanks.

Answer 1

In step $5$ you are going to multiply the digit $d$ you put on top by $20+d$ because you write $d$ to the right of the $2$ and multiply. To get a positive remainder you need $d(20+d) \le 174$. Here if $d=6, d(20+d)=156$, while if $d=7, d(20+d)=189,$ which is too large.

The idea of the algorithm is based on $(a+b)^2=a^2+2ab+b^2$. $a$ represents the digits you have already found, and you have subtracted $a^2$ from the number you are taking the square root of already. In your example, the first $a$ is $100$ and we subtracted $100^2$ in step $3$. Writing twice the current set of digits (here $2$) is the $2a$ part of $2ab+b^2$. When you write the next digit in both places and multiply, you get $2ab+b^2$. When you subtract that you have subtracted $a^2+2ab+b^2$ from the original number, here leaving $27493-25600=1893$. Our new $a$ is $160$ and we need to find $b$ such that $2ab+b^2 \lt 1893$, which turns out to be $5$.

Answer 2

In case you read algorithms better than words (nobody has to be that geeky :), here is an algorithm in pseudo-code that does what is described on Wikipedia (from the comment by @JohnOmielan). Unfortunately, Wikipedia doesn't offer an algo.

  input ( n ) 
  mantissa = mantissa of n
  decimals = decimals of n

  # prefix a "0" to mantissa or postfix "0" to decimals if necessary
  if (length (mantissa) is odd) mantissa = "0" + mantissa
  if (length (decimals) is odd) decimals = decimals + "0"

  # rebuild n                                                 
  n = mantissa + "." + decimals

  # empty result, remainder, and p (root) variables
  result = ""
  remain = 0
  p = 0

  # loop 
  while ( n has characters )
  do
     a = takeFirstCharacter( n )
     if ( a = "." )
      then
          addDecimalPointTo( result )
      else
          b = takeFirstCharacter( n )
          c = remain*100 + a*10 + b
          # find smaller x such that x*(20*p + x) <= c
          x = 0
          while ((x+1)*(20*p + (x+1)) <= c)
          do
             x = x + 1
          done
          y = x*(20*p + x)
          p = 10*p + x      # "push" x at the end of p
          remain = c - y    # remainder
          concatenateDigit (x) to ending of (result)
      end if
  end while

  # 'result' is the sqrt

For those interested, below is an actual C program that implements the above algorithm.

An actual program should be quite straightforward, but taking into account the special cases (odd lengths, no decimals/mantissa...) and the required precision, the code is not so short ; the effective loop itself is short though.

This is an $O(length(n))$ algorithm (complexity driven by the Wikipedia article).

Of course, the intermediary calculations (like the current root that grows by a factor $10$ at each iteration, even after the decimal dot) requires to select the bigger integer type as possible (I chose (gcc) $T=long \ long$, or $64$ bits).

Also, the temporary strings, $n$ with padded $0$es for oddies and precision, and result, must have sufficient lengths to bear the numbers provided.

Result is a string (not a double), it's simpler to deal with, and, at least, any precision loss will not be due to the floating point conversions and operations (...).

Can be optimized (the multiple strcat for instance).

First, a type definition and the conversion macros

typedef long long T; // type to use for growing integers

// Conversion from char to digit, and reverse
#define C2D(c)      ((c)-'0')
#define D2C(d)      ((d)+'0')

The function that takes the $arg$ argument and precision

char *sqr(char *arg, int precision) {
     static char res[100];   // assuming result (watch precision) won't exceed 99 chars
     char n[100] = { 0 };    // assuming input after adjustments won't exceed 99 chars

     // Calc lengths of mantissa and dec parts
     int len = strlen(arg);
     char *dot = strchr(arg, '.');                  // do we have a '.'?
     int dotpos = dot ? (int)(dot - arg) : len;     // '.' position (or end of n)
     int declen = dot ? len - dotpos - 1 : 0;       // decimals part length

     // Mantissa
     if ( ! dotpos) strcat (n, "00");       // no mantissa, set it at "00"
     else if (dotpos & 1) strcat (n, "0");  // mantissa odd length, prefix a "0"
     dotpos += strlen(n);   // adjust dot position if "0"oes were added
     
     strcat(n, arg); // concat the argument to n (being "", "0" or "00")
     
     // Decimals
     if (declen > precision*2) n[dotpos+1+precision*2] = 0; // truncate decs?
     else {
          if (dot && (declen & 1)) {
                strcat(n, "0"); // decimals odd length, postfix "0"
                declen++;
          }
          // postfix as many "00" as necessary to reach "precision"
          if (declen < precision*2) {
                if ( ! dot) strcat(n, "."); // ensure we have a '.'
                while (declen < precision*2) {
                     strcat(n, "00");
                     declen += 2;
                }
          }
     }

     // Init algo
     char *s = n;           // points to first char
     char *r = res;         // points to result
     T remain = 0;          // remainder
     T p = 0;               // root so far

     // Loop until no more chars
     while ( *s ) {
          if (*s == '.') {  // found the '.'
                *r++ = '.'; // store a '.' in result
                s++;        // advance our string cursor by 1
          }
          else {
                T c = remain*100 + C2D(*s)*10 + C2D(s[1]); // remain*100 + 2-digits
                // find smaller x for that equation to be <= c
                int x = 0;
                while((x+1)*(20*p + (x+1)) <= c) x++;
                T y = x*(20*p + x);
                p = 10*p + x;           // our root
                remain = c - y;         // new remainder
                *r++ = D2C(x);          // store next root digit (char)
                s += 2;                 // advance our string position in n
          }
     }
     *r++ = 0; // "close" the string (in C, "strings" end with a byte 0)

     return res;
}

That might be called from main() this way

int main(int argc, char **argv){

     // Usage: ./prog n [precision]

     char *arg = *++argv;
     int precision = 3; // 3 dec precision by default
     if (argc > 2) precision = atoi(*++argv);

     printf("sqrt(%s) = %s\n", arg, sqr(arg, precision));

     return 0;
}