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 <float.h>

#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 <math.h>
#include <errno.h>

#include "decompose.h"

TYPE NAME_SFX(sqrt)(TYPE a)
 (r.o.e >> 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.

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*)&f; // double to int
i = i - (1023 << 52); // remove IEEE754 bias
i = i >> 1; // divide exponent by 2
i = i + (1023 << 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*)&f; // double to int
i = i >> 1; // divide exponent by 2
i = i + (1023 << 52 - 1023 << 51); // remove/add IEEE754 bias

The constant can then be simplified to:

uint64_t i = *(uint64_t*)&f; // double to int
i = (i >> 1) + (1023 << 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. 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 mapping 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 )

My mapping is far more accurate and yours isn't even monotonic.

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*)&x;
 *i = (*i >> 1) + (127 << 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*)&x;
 *i = (*i >> 1) + (1023 << 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;

Modular odd check

The compiler is likely to do this anyway, but can't

if(c % 2 == 0)

be

if (!(c & 1))

?

Factoring

Isn't

r.f -= (r.f*r.f-a)/(2*r.f);

equivalent to:

r.f = (r.f + a/r.f)/2?

Optimization

Depending on which compiler you're using and what flags you pass it, some of your attempted math risks modification by the optimizer, and in an implementation so closely coupled to the hardware, this might matter. Have you checked your produced assembly? Can you show us your build flags?

I don't see why you define this constant yourself:

#define MANTISSA_SIZE 52

Given we already assume that FLT_RADIX is 2, we can use the appropriate macro from <float.h> (DBL_MANT_DIG for double, etc.).

I think there's danger of integer overflow here:

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

We'd be better extracting into a big enough temporary, and applying the exponent bias to that:

int e = r.o.e - EXPONENT_BIAS;
e = (e & 1) | (e >> 1);
r.o.e = e + EXPONENT_BIAS;

It might be possible to shift in place and then correct using half the bias; I haven't checked whether that's equivalent.