sqrt.c File Reference

#include <math.h>
#include <assert.h>

Include dependency graph for sqrt.c:

Go to the source code of this file.

Functions
double __cdecl	sqrt (double x)

Function Documentation

sqrt()

double __cdecl sqrt

(

double

x

)

Definition at line 13 of file sqrt.c.

{
    const double threehalfs = 1.5;
    const double x2 = x * 0.5;
    long long bits;
    double inv, y;
 
    /* Handle special cases */
    if (x == 0.0)
    {
        return x;
    }
    else if (x < 0.0)
    {
        return -NAN;
    }
 
    /* Convert into a 64  bit integer */
    bits = *(long long *)&x;
 
    /* Check for !finite(x) */
    if ((bits & 0x7ff7ffffffffffffLL) == 0x7ff0000000000000LL)
    {
        return x;
    }
 
    /* Step 1: quick approximation of 1/sqrt(x) with bit magic
       See http://en.wikipedia.org/wiki/Fast_inverse_square_root */
    bits = 0x5fe6eb50c7b537a9ll - (bits >> 1);
    inv = *(double*)&bits;
 
    /* Step 2: 3 Newton iterations to approximate 1 / sqrt(x) */
    inv = inv * (threehalfs - (x2 * inv * inv));
    inv = inv * (threehalfs - (x2 * inv * inv));
    inv = inv * (threehalfs - (x2 * inv * inv));
 
    /* Step 3: 1 additional Heron iteration has shown to maximize the precision.
       Normally the formula would be: y = (y + (x / y)) * 0.5;
       Instead we use the inverse sqrt directly */
    y = ((1 / inv) + (x * inv)) * 0.5;
 
    //assert(y == (double)((y + (x / y)) * 0.5));
    /* GCC BUG: While the C-Standard requires that an explicit cast to
       double converts the result of a computation to the appropriate
       64 bit value, our GCC ignores this and uses an 80 bit FPU register
       in an intermediate value, so we need to make sure it is stored in
       a memory location before comparison */
//#if DBG
//    {
//        volatile double y1 = y, y2;
//        y2 = (y + (x / y)) * 0.5;
//        assert(y1 == y2);
//    }
//#endif
 
    return y;
}