static unsigned long isqrt(const unsigned long X)
{
    unsigned long k = 1, result = 1;
    while (result <= X)
    {
        k++;
        result = k*k;
    }
    return k--;
}

static double nsqrt(double X)
{
    const double inf = 1.0l/0.0l; /* positive infinity */
    const double ninf = -inf; /* negative infinity */
    const double nan = 0.0l/0.0l; /* not a number */
    unsigned long I;
    double A_prev, A;
    if (X == 0.0 || X == 1.0 || X == inf)
    {
        return X;
    }
    else if (X == 2.0)
    {
        return M_SQRT2;
    }
    else if (X == ninf || X < -0.0)
    {
        return nan;
    }
    /* NaN cannot be compared using == operator */
    else if (memcmp((void *)&X, (void *)&nan, sizeof(double)) == 0)
    {
        return nan;
    }
    /* I =  floor(sqrt(X)) */
    I = isqrt((unsigned long)X);
    if (X == (double)(I*I)) /* exact root */
    {
        return (double)I;
    }
    /* Newton's method: A{n} = (A{n-1} + X/A{n-1})/2 */
    A = A_prev = I; /* A{1} = I, or A{1} = floor(sqrt(X)) */
    do {
        A_prev = A;
        A = (A + X/A)/2.0;
    } while (A_prev != A);
    return A;
}