# Approximating the square root using an iterative method

I wrote this code, based on the Newton-Raphson method, to find the square root of a number. I'm wondering how I can optimise this code, as I am out of ideas.

#include <stdio.h>

int main(void)
{
float n, x, i;

printf("Enter the number you wish to find the square root of.\n");
printf("\n");
scanf("%f", &n);

x = n/2;

for (i = 0; i < 100; i++)
x = x - (((x*x) - n)/(2*x));

printf("The square root of %.0f is %.4f.\n", n, x);

return 0;

}

• Better to use %e for output. Else small values will all have a square root of 0.0000 and large number will print with superfluous digits. – chux May 2 '14 at 13:23
• Or just include math.h and use sqrt... What do you say? – user3459110 May 2 '14 at 15:07
• The thing is, I'm required to write a function that finds an approximation to the square root of a number, without using the sqrt function. – user41697 May 3 '14 at 9:32

• Do not rely on a fixed number of iterations. Stop when reaching a desired discrepancy (e.g. fabs(x*x - n) < epsilon)

• I do not see any reason for not simplifying the calculations into x = (x + n/x)/2.0.

• Assure that the input is positive.

• Thanks for this! Just one question, what's the value of epsilon or how do I choose it? – user41697 May 2 '14 at 0:09
• @user3182162: With floating-point numbers, it's based on the size of the number. Assuming n is positive, a useful epsilon would be somewhere around pow(2.0, floor(log2(n)) - 23). (Of course, that might not help much if you're trying to find the square root iteratively -- the square root itself is pow(n, 0.5), so the iteration almost wouldn't be worth the trouble. :P) – cHao May 2 '14 at 0:54
• Floating point numbers have a logarithmic distribution. Stopping on a delta value like epsilon is a poor approach. Consider epsilon = 1e-6 and performing square_root(1e20) or square_root(1e-20). Better to use a relative__epsilon e.g. fabs(x*x - n) < (epsilon*n). Even better, given the convergence rate of this algorithm, I believe it converges in about ln2(significand_bits_in_float) (e.g. 5 or 6) iterations. – chux May 2 '14 at 13:08

The square root procedure could easily be extracted into its own function, and therefore should be.

i is a counter, and therefore should be an int, not a float. As @fscore points out, 100 iterations is overkill. Assuming that each iteration gains you an extra bit of precision — in other words, halves the error — there aren't 100 bits in a float! And if you really cared about precision, you should be using doubles everywhere.

• +1 for i. This method doubles the number of corrects bit per iteration. So with float about 24 bits, ln2(24) --> 5. So about 5 iterations needed. – chux May 2 '14 at 13:15

As you probably know, once it iss close, Newton-Raphson converges really quickly.

Therefore, if you start from a reasonable approximation, you need a limited number of iterations.

If you represent n as as power of 2 times a fraction 1.?????, then the square root will be given by half that power of 2 times a fraction.

So, if you can find that power, you have a good start.

Now, that is exactly how a floating point number is represented internally. See Wikipedia for the IEEE 754 double-precision binary floating-point format.

We can therefore write:

#include <stdint.h>

const int signPosition = 63;
const int fractionBits = 52;

const uint64_t signMask = 1ull << signPosition;
const uint64_t fractionMask = (1ull << fractionBits) - 1;
const int exponentBias = 1023;

int getExponent(double x)
{
uint64_t representation = *(uint64_t*)&x;
uint64_t exponentBiased = (representation & exponentMask) >> fractionBits;
return (int)(exponentBiased)-exponentBias;
}

double firstGuessOfSquareRoot(double x)
{
unsigned exponent = getExponent(x);
unsigned rootExponent = exponent / 2;
double guess = 1 << rootExponent;
return guess;
}

double root(double x)
{
double guess = firstGuessOfSquareRoot(x);
int loops = 6;
while (loops > 0)
{
guess += (x - guess * guess) / 2 / guess;
--loops;
}
return guess;
}


I'd suggest 6 loops is not too bad.

Now, you can do even better, getting down to 4 or maybe 3 loops if you start with a better guess. That can be done by:

• Looking at the fractional part of the double and using a look-up table.
• Handling the case where the original exponent is odd correctly - the above gives an initial guess out by at least sqrt(2) for the case.

Oh - some might ask: Why use a fixed number of iterations? The reason is that you can show that no more than this number is required and also that this number is typically required. Adding an error margin based check will slow down the main case whilst certainly being faster for special cases. I'm not sure as to the tradeoff.

• +1 for good convergence discussion and better for fix loop count (integer math) rather than FP termination calculation. – chux May 2 '14 at 13:17
• +1 for convergence discussion. For production code, I'd prefer defining getExponent() using the standard C library frexp() function rather than re-writing yet another redundant implementation of it. – David Cary May 2 '14 at 15:38
• @DavidCary. Didn't know about frexp() - agree one should use if possible. But then, we are intentionally reinventing wheels here, so I'll leave the code as is. – Keith May 4 '14 at 23:41
• Additional newlines can just be in the first output statement. You should also use puts() instead since this is unformatted output. It also adds an additional newline.

puts("Enter the number you wish to find the square root of.\n");

• I'd recommend renaming n and x to something like input and result respectively.

Avoid using single-character variables unless they're for loop counters.

• First point: the compiler will look after that automatically. – edmz May 2 '14 at 8:54

There are a couple of possible ways to optimize the code. @Keith has already done what I'd consider quite a good job of covering the micro-optimizations within this algorithm by improving your initial guess and your loop exit criteria.

The other obvious possibility would be to use a different algorithm. For one example, the binary reducing algorithm can compute square roots to the nearest integer quite a bit faster than the code you posted, and still somewhat faster than the code @keith posted.

int isqrt(int input) {
int result = 0;
int bit = 1 << (sizeof(int)*CHAR_BIT-2);

while (input <= bit)
bit /= 4;

do {
if (input >= result + bit) {
input -= result + bit;
result = result / 2 + bit;
}
else
result /= 2;
bit /= 4;
} while (bit != 0);

return result;
}


The primary advantage of the binary reducing method is that the divisions it uses are both by fixed powers of two (specifically, divisions by 2 and 4). On a binary computer with a halfway decent compiler, you can pretty well count on these being implemented as bit shifts.

Translation to floating point isn't quite so straightforward though. In a Newton-Raphson implementation, a better initial guess improves speed, but (at least for square root) basically can't ever result in incorrect results. By contrast, the:

int bit = 1 << (sizeof(int)*CHAR_BIT-2);

while (input <= bit)
bit /= 4;


...part of this algorithm must be correct to produce correct results at all. Specifically, it needs to be the largest power of 4 that's less than or equal to the input. You obviously can compute that value in floating point, but doing so is somewhat non-trivial compared to doing it for an integer.

Newton-Raphson will usually require fewer iterations, but each iteration is slower because its divisions are by numbers generated from the input, so they need to be executed as real division operations. For integers this gains speed because although it needs more iterations, each iteration is quite fast. In floating point, you end up using a real divide instruction for each of those division operations, losing most of the potential speed gain. If you're willing to do like @Keith did and directly manipulate the bits of the floating point representation, you might be able to gain this back, because the divisions can actually be simple subtractions from the exponent part of the floating point number.

Bottom line: this produces good results fairly quickly for integers. Conversion to floating point is possible, but somewhat non-trivial. Doing so in a way that retains much (if any) speed advantage starts to move into the decidedly non-trivial range.

• Whoa! Thanks dude, quite a detailed description, and easy to follow concepts. The kind of answers that really help! – user41697 May 4 '14 at 17:18

Some recommendations and changes in code,

• better variable names
• Prompting user if he wants to continue or stop the program preventing infinite loop
• Formula reduced to simple easy to understand format

#include<stdio.h>

int main(void)
{
float input_num, output; //better variable names
int continue_input=1;

while(continue_input == 1) //user input continuous
{
puts("Enter the number you wish to find the square root of.\n");
scanf("%f", &input_num);

output = input_num/2;

output *= (2 - input_num)/2;

printf("The square root of %.0f is %.4f.\n", input_num, output);

puts("Enter 0 to stop & 1 to continue: ");
scanf("%d",&continue_input);
}
return 0;
}


Why do you need 100 iterations? They seem really unnecessary as you are taking user input and i is not used anywhere, but again you need to check that out.

• You can use output *= ... instead of output = output * .... – Jamal May 2 '14 at 2:25
• This also appears to be an infinite loop (specifically while (1)). val is never changed within the loop. You could ask the user if they'd like to continue, and break if it's a no. If you go with that, then you can just have while (1) and remove val (which isn't a very good name anyway). – Jamal May 2 '14 at 2:34
• @Jamal check the changes – fscore May 2 '14 at 2:46
• Looks better. You can also declare i right before the for loop to keep it in close scope (assuming the OP is not using C99). – Jamal May 2 '14 at 2:52
• The second line of calculations with output still doesn't work, based on my tests on Ideone.com. The OP's does work, but I know you're trying to simplify it. – Jamal May 2 '14 at 3:10