IEEE 754 square root with Newton-Raphson I have an implementation of the sqrt function that uses a combination of IEEE 754, packed bitfields, and the Newton-Raphson algorithm:

decompose.h:

#ifndef DECOMPOSE_H
#define DECOMPOSE_H 1

#ifdef __cplusplus
extern "C" {
#endif

#include &lt;float.h&gt;

#define MANTISSA_SIZE 52
#define EXPONENT_SIZE 11
#define EXPONENT_BIAS 0x3ff
#define TYPE double
#define TYPE_MIN DBL_MIN
#define TYPE_MAX DBL_MAX
#define NAME_SFX(name) name

typedef unsigned long long mantissa;
typedef unsigned short exponent;
typedef unsigned char sign;

#include "decompose_generic.h"

#ifdef __cplusplus
}
#endif

#endif


decompose_generic.h:

#ifndef DECOMPOSE_GENERIC_H
#define DECOMPOSE_GENERIC_H 1

#ifdef __cplusplus
extern "C" {
#endif

#define SIGN_SIZE 1

union decompose {
    TYPE f;
    struct decompose_s {
        mantissa m:MANTISSA_SIZE;
        exponent e:EXPONENT_SIZE;
        sign s:SIGN_SIZE;
    } __attribute__((packed)) o;
};

static inline union decompose decompose(TYPE f)
{
    union decompose result;
    result.f = f;
    return result;
}

#define is_nan(x) NAME_SFX(isnan)(x)

#ifdef __cplusplus
}
#endif

#endif


sqrt.c:

#include &lt;math.h&gt;
#include &lt;errno.h&gt;

#include "decompose.h"

TYPE NAME_SFX(sqrt)(TYPE a)
{
    TYPE result_old2 = 0;
    unsigned short c = 0;
    union decompose r = decompose(a);

    if(a &lt; -(TYPE)0)
    {
        errno = EDOM;
        return NAME_SFX(nan)("");
    }
    if(a == (TYPE)0 || is_nan(a) || a &gt; TYPE_MAX)
    {
        return a;
    }

    /* Divide the exponent by 2 */
    r.o.e -= EXPONENT_BIAS;
    r.o.e = (r.o.e &amp; 1) | (r.o.e &gt;&gt; 1);
    r.o.e += EXPONENT_BIAS;

    /* Newton-Raphson */
    while(result_old2 != r.f)
    {
        if(c % 2 == 0) result_old2 = r.f;
        r.f -= (r.f*r.f-a)/(2*r.f);
        c++;
    }

    return r.f;
}


What I'm looking for:


Are there any efficiency problems?
Are there any accuracy problems?
Are there any other miscellaneous things that I can improve?


As far as I've tested, this is as accurate as possible. It even matches the results from the GNU C library.

One note: This is part of a modular system. There will be a decompose.h for each floating point type (float, double, long double), and they all include decompose_generic.h.

I'm slightly new to this area of C (math, IEEE 754), but not so new that I'd add the beginner tag.
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.