# Help to understand algorithm for finding square root?

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.

• FYI, the algorithm is described (IMHO, somewhat better in at least certain ways, including the Step #$5$ as explained in Ross' answer below, with it being the second bullet point of its step #$2$) in the Decimal (base 10) section of Wikipedia's "Methods of computing square roots" article. Jun 16, 2020 at 4:26
• Also, FYI, the question at Derivation of Square-Root Computation Algorithm asks to have explained the Wikipedia algorithm section I link to in my comment above. Jun 16, 2020 at 4:50

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$$.

• Thanks, that was very clarifying =)
– PtF
Jun 16, 2020 at 12:56
• Now I see the confusion, I was taking "the smallest remainder" as the remainder of a division, maybe it is a language issue since english is not my mother tongue.
– PtF
Jun 16, 2020 at 12:59
• I think the language is confusing. You want to choose $d$ so you get the smallest nonnegative result from the subtraction on that line. They are calling the result the remainder. Jun 16, 2020 at 13:39

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
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;   // assuming result (watch precision) won't exceed 99 chars
char n = { 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"

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); // 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;
}

• Algorithms should be presented in pseudo-code... This will be helpful to more people. Also, pseudo-code clears the syntax noise and language-dependent subtleties away. Jun 16, 2020 at 8:55
• The C program was intended for those who want to play with the code, but, indeed, the pseudo-code is more meaningful in this case, and on Maths SE in general. It has been added near the top. Jun 16, 2020 at 9:33
• Thanks, I'll try to implement it =)
– PtF
Jun 16, 2020 at 12:55