FYI, in IEEE754 sqrt is a "basic" operation that's required to be correctly-rounded (rounding error <= 0.5ulp), same as + - * /. Hardware FPUs (I think) always provide sqrt if they provide the other operations, especially division which is typically implemented similarly (and with similar performance).
NR iterations each involving a division are not going to be faster than HW sqrt. I think that's part of why x86's approximation instructions are for reciprocal and rsqrt, not sqrt: the Newton iteration formula for 1/sqrt(x) doesn't involve division, so on old CPUs with very slow sqrtps compared to mulps, it was worth using sqrt(x) ~= x * approx_rsqrt(x), optionally with a Newton iteration on the approx_rsqrt result. Maybe sometimes still worth doing that in a loop that does nothing else, but normally you should aim for more computational intensity (do more with your data while it's hot in registers; don't write loops that just load/sqrt/store).
IDK if it's possible for your code to converge on a value that's not the correctly-rounded result, but something you should test for by checking against the implementation's sqrt function.
(Call your function something else. With GCC for example, sqrt is a special name unless you use -fno-builtin-sqrt. Callers will inline the sqrt instruction on targets that have one, like x86 sqrtsd, only calling the libc function for inputs that need to set errno. Or never if you compile with -fno-math-errno, always a good idea. Math-errno is an optional feature that not all libcs support, although GNU C does. It's a pretty bad obsolete design that nobody should use; that's why we have fenv.h to check per-thread sticky FP exception flags.)
Wrong initial guess for inputs < 1.0
/* 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;
Your exponent bitfield type is unsigned short1 so this is a logical right shift.
But anyway, logical right shift is a bug in your initial approximation.
The exponent field is a 2's complement integer (encoded with a bias, but that's just an encoding/decoding detail you already take care of). A logical right shift of a small exponent (like for 0.0001) will result in a large positive exponent from shifting in a zero. Your initial guess for 0.0001 will be something close-ish to DBL_MAX, way above 1.0.
You need an arithmetic right shift of the exponent to bring it closer to the middle exponent value (unbiased) 0 or (biased) 127. Thus bringing the magnitude of the represented value closer to 1.0.
Most C implementations make the sane choice that >> of a signed integer type is an arithmetic right shift. ISO C requires that implementations define which it is, but unfortunately leaves them the choice of logical or arithmetic. (Unlike signed left shift, it can't be UB).
But fortunately your code is written in GNU C (not ISO C), e.g. using __attribute__((packed)) on a struct. GNU C guarantees that >> on signed integers is arithmetic. (And that signed integer types are 2's complement!) So the only implementations that can compile your code are ones where >> on an int does what you want.
(GCC manual, integer implementation-defined behaviour):
Signed ‘>>’ acts on negative numbers by sign extension.
You might want to use unsigned bitfield members and unsigned arithmetic for bias/unbias. But you can safely use int32_t (which is guaranteed to be a 2's complement integer type), or just int for the temporary holding the field value for the bias/unbias calculations. Bias/unbias may cross zero (unsigned wrapping), not signed overflow.
Or better, use @Rainer P's method:
Don't unbias first, just unsigned shift and then add half the bias constant. (See his answer for how it simplifies down to that). Doing the shift before unbiasing means the exponent is unsigned so shifting in a zero is good. (I think this works right even for tiny values with biased-exponents near zero, i.e. near the smallest possible exponent encoding).
I tried a few values on https://www.h-schmidt.net/FloatConverter/IEEE754.html, manually shifting the exponent and then adding 1 (with carry) to the 2nd-highest exponent bit. e.g. 2.3509887E-38 (biased exponent = 0b0000'0010) becomes 2.1684043E-19 (biased exponent = 0b0100'0001).
But beware when shifting the entire bit-pattern that you don't shift in a 1 from the sign bit of -0.0. You are handling a == 0.0 as a special case so you're fine there because -0.0 == 0.0 is true. I guess you can't let 0.0 go through the main NR loop because you'd have a division by zero.
Footnote 1:
IDK why you'd choose unsigned short instead of just unsigned for your exponent bitfield type. The underlying type of a bitfield can be wider than the field without hurting anything, and some compilers will do math at the width of that type. So using unsigned int if it's wide enough is a good idea. int is at least 16 bits.
Using uint64_t as the bitfield type can result in slow 32-bit code with some compilers (notably MSVC), but I think unsigned int is always fine for the exponent. Unsigned does make sense for the biased exponent.
Fun fact: just int as the bitfield type leaves the signedness up to the implementation. You can use signed int to force signed. https://stackoverflow.com/questions/42527387/signed-bit-field-represetation
unroll instead of checking c % 2 == 0
A compiler might do this anyway, but it probably makes more sense to just unroll your loop by 2 instead of using a counter to do something special every other iteration.
The actual Newton iteration expression is compact enough that repeating it once is probably more simple for humans to read than working through the c++ and if (c%2 == 0) logic.
You can handle the while(foo) condition as if(!foo) break; when unrolling, or simply leave it out and only check for leaving the loop every 2 Newton iterations, unless you need to check both possibilities to catch alternation between two values.
(With a bad initial estimate that will take many iterations anyway, this might be worth it. Otherwise probably not, especially on a superscalar / out-of-order CPU that can be checking the branch condition in parallel with the main dependency chain which doesn't have a lot of ILP (instruction-level parallelism); mostly one long dependency chain.)
Pull some checks off the fast-path
Your version has 4 separate checks on a before we reach the actual work.
if(a < -(TYPE)0)
{
errno = EDOM;
return NAME_SFX(nan)("");
}
if(a == (TYPE)0 || is_nan(a) || a > TYPE_MAX)
{
return a;
}
Instead, catch the special cases with one or two compares, and sort them out in that block, off the fast path. The performance of sqrt(NaN) is much much less important than the performance of sqrt(1.234)
if(! a > (TYPE)0)
{
if(a == (TYPE)0 || is_nan(a))
return a;
// else a < 0
errno = EDOM;
return NAME_SFX(nan)("");
}
if( a > TYPE_MAX)
{
return a;
}
If you need to handle subnormals specially (at least for generating the initial approximation?), you might use if (! a>=DBL_MIN). Or you could decide that your soft-float implementation simply doesn't handle subnormals.
IDK what your goal is with this code; if it was part of a soft-float library, simply using FP add, mul, and divide seems inefficient. You'd want to optimize multiple operations together and so on.
But if it's for HW FPUs, then I think real all FPUs will have sqrt built-in as well.
So I assume this is just a learning exercise. So we won't get into too detail of optimizing it for any particular hardware or compiler.
std::numeric_limits<T>::is_iec559. I don't know of an equivalent test for C, other than separately checking all the macros in<float.h>. \$\endgroup\$__STDC_IEC_559__\$\endgroup\$