Methods of Computing Square Roots - Approximations That Depend On IEEE Representation

Approximations That Depend On IEEE Representation

On computers, a very rapid Newton's-method-based approximation to the square root can be obtained for floating point numbers when computers use an IEEE (or sufficiently similar) representation.

This technique is based on the fact that the IEEE floating point format approximates base-2 logarithm. For example, you can get the approximate logarithm of 32-bit single precision floating point number by translating its binary representation as an integer, scaling it by, and removing a bias of 127, i.e.

For example, 1.0 is represented by a hexadecimal number 0x3F800000, which would represent if taken as an integer. Using the formula above you get, as expected from . In a similar fashion you get 0.5 from 1.5 (0x3FC00000).

To get the square root, divide the logarithm by 2 and convert the value back. The following program demonstrates the idea. Note that the exponent's lowest bit is intentionally allowed to propagate into the mantissa. One way to justify the steps in this program is to assume is the exponent bias and is the number of explicitly stored bits in the mantissa and then show that

/* * Warning: the following program has undefined behavior. * Reading u.tmp after storing u.f is not allowed by the C standard, nor the C++ standard. */ float sqrt_approx(float z) { union { int tmp; float f; } u; u.f = z; /* * To justify the following code, prove that * * ((((val_int / 2^m) - b) / 2) + b) * 2^m = ((val_int - 2^m) / 2) + ((b + 1) / 2) * 2^m) * * where * * val_int = u.tmp * b = exponent bias * m = number of mantissa bits * * . */ u.tmp -= 1 << 23; /* Subtract 2^m. */ u.tmp >>= 1; /* Divide by 2. */ u.tmp += 1 << 29; /* Add ((b + 1) / 2) * 2^m. */ return u.f; }

The three mathematical operations forming the core of the above function can be expressed in a single line. An additional adjustment can be added to reduce the maximum relative error. So, the three operations, not including the cast, can be rewritten as

u.tmp = (1 << 29) + (u.tmp >> 1) - (1 << 22) + a;

where a is a bias for adjusting the approximation errors. For example, with a = 0 the results are accurate for even powers of 2 (e.g., 1.0), but for other numbers the results will be slightly too big (e.g.,1.5 for 2.0 instead of 1.414... with 6% error). With a = -0x4C000, the errors are between about -3.5% and 3.5%.

If the approximation is to be used for an initial guess for Newton's method to the equation, then the reciprocal form shown in the following section is preferred.