There is a better algorithm, which needs at most 6 iterations to converge to maximum precision for double numbers:
#include <math.h> double sqrt(double x) { if (x <= 0) return 0; // if negative number throw an exception? int exp = 0; x = frexp(x, &exp); // extract binary exponent from x if (exp & 1) { // we want exponent to be even exp--; x *= 2; } double y = (1+x)/2; // first approximation double z = 0; while (y != z) { // yes, we CAN compare doubles here! z = y; y = (y + x/y) / 2; } return ldexp(y, exp/2); // multiply answer by 2^(exp/2) }
Algorithm starts with 1 as first approximation for square root value. Then, on each step, it improves next approximation by taking average between current value y
and x/y
. If y
= sqrt(x)
, it will be the same. If y
> sqrt(x)
, then x/y
< sqrt(x)
by about the same amount. In other words, it will converge very fast.
UPDATE: To speed up convergence on very large or very small numbers, changed sqrt()
function to extract binary exponent and compute square root from number in [1, 4)
range. It now needs frexp()
from <math.h>
to get binary exponent, but it is possible to get this exponent by extracting bits from IEEE-754 number format without using frexp()
.