IEEE 754 square root with Newton-Raphson
This answer uses pointer-casting for type-punning just to save space.  In practice keep using your union (safe in ISO C99, and in C++ as a GNU and MSVC extension) or memcpy (safe in C and C++).  This pointer-casting is only safe in MSVC, or GNU-compatible compilers with -fno-strict-aliasing

Initial approximation

Packed bit fields are not only unnecessary here, they make the result worse. You do something along the lines of:

uint64_t i = *(uint64_t*)&amp;f;           // double to int
i = i - (1023 &lt;&lt; 52);                  // remove IEEE754 bias
i = i &gt;&gt; 1;                            // divide exponent by 2
i = i + (1023 &lt;&lt; 52);                  // add IEEE754 bias


If you change the order of operations, you can handle the biasing in a single constant:

uint64_t i = *(uint64_t*)&amp;f;           // double to int
i = i &gt;&gt; 1;                            // divide exponent by 2
i = i + (1023 &lt;&lt; 52) - (1023 &lt;&lt; 51);   // remove/add IEEE754 bias


The constant can then be simplified to:

uint64_t i = *(uint64_t*)&amp;f;           // double to int
i = (i &gt;&gt; 1) + (1023 &lt;&lt; 51);           // all of the above


Note that I didn't bother to mask when shifting, so the exponent's rightmost bit drops into the mantissa, which is shifted as well. This is a feature, not a bug. The IEEE754 formats are deliberately chosen to be monotonic when interpreted as integers, so carry-overs between exponent and mantissa are explicitly allowed.

With my approximation, the mapping is as follows (linear mapping within intervals):

original   number is    exponent   manitssa   orig interval
exponent  in interval   in sqrt   first bit     maps to
   0     [  1.0,  2.0 )    0         0       [  1.0,  1.5 )
   1     [  2.0,  4.0 )    0         1       [  1.5,  2.0 )
   2     [  4.0,  8.0 )    1         0       [  2.0,  3.0 )
   3     [  8.0, 16.0 )    1         1       [  3.0,  4.0 )
   4     [ 16.0, 32.0 )    2         0       [  4.0,  6.0 )


With your code, the mapping is (also linear within intervals):

original   number is    exponent  orig interval
exponent  in interval   in sqrt     maps to
   0     [  1.0,  2.0 )    0     [  1.0,  2.0 )
   1     [  2.0,  4.0 )    0     [  1.0,  2.0 )
   2     [  4.0,  8.0 )    1     [  2.0,  4.0 )
   3     [  8.0, 16.0 )    1     [  2.0,  4.0 )
   4     [ 16.0, 32.0 )    2     [  4.0,  8.0 )


Newton-Raphson

Given a good approximation, Newton-Raphson doubles the number of significant digits on each iteration (quadratic convergence). The above approximation provides about 4 bits of accuracy (max error: 6% or ~1/16), so 3 Newton-Raphson iterations are required for single and 4 iterations for double precision.

The entire code could look like this:

float sqrtf(float x)
{
    float y = x;
    // Approximation
    uint32_t* i = (uint32_t*)&amp;x;
    *i = (*i &gt;&gt; 1) + (127 &lt;&lt; 22);
    // Newton-Raphson
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    return x;
}

double sqrt(double x)
{
    double y = x;
    // Approximation
    uint64_t* i = (uint64_t*)&amp;x;
    *i = (*i &gt;&gt; 1) + (1023 &lt;&lt; 51);
    // Newton-Raphson
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    x = (x + y/x) / 2;
    return x;
}


Subnormals and other hazards

You already have special handling for nan, inf, 0 and x &lt; 0, which is good, but we need to handle subnormal numbers as well. The code, both yours and mine, fails when confronted with such numbers. As a remedy, I suggest you scale subnormal numbers by 2^200 and scale the result by 2^-100. I chose 2^200 because it safely brings all subnormals back into normal range for all precisions, including quad, but other numbers work equally well (but only powers of two retain full precision).

Also observe that my code relies on IEEE754 properties. It doesn't work for non-IEEE754 numbers, most notably x86's 80bit extended double format. If you want to handle such numbers, you can convert them to float or double for the approximation step and back to extended double for Newton-Raphson.