23
\$\begingroup\$

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)
{
    TYPE result_old2 = 0;
    unsigned short c = 0;
    union decompose r = decompose(a);

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

    /* 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;

    /* 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 tag.

\$\endgroup\$
  • 9
    \$\begingroup\$ It's a shame this isn't C++, where you could check that the floating-point representation matches your expectation using 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\$ – Toby Speight Sep 24 at 16:39
17
\$\begingroup\$

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*)&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. 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*)&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;
}

Subnormals and other hazards

You already have special handling for nan, inf, 0 and x < 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.

\$\endgroup\$
  • \$\begingroup\$ This won't work for 80-bit long double, which can't be represented as an integer. \$\endgroup\$ – JL2210 Sep 24 at 17:59
  • 2
    \$\begingroup\$ Also note that this violates strict aliasing and is undefined behavior. \$\endgroup\$ – JL2210 Sep 24 at 19:39
  • 4
    \$\begingroup\$ Your approximation is off by 41% at worst, that's about 1bit. (with my approximation it's ~6%) - @JL2210 \$\endgroup\$ – Rainer P. Sep 24 at 19:48
  • 5
    \$\begingroup\$ *(uint64_t*)&f is undefined behaviour in ISO C. It's only well-defined in MSVC or gcc -fno-strict-aliasing. , @JL2210 used a union to type-pun because that's well-defined in ISO C99. (But in C++ only as a GNU extension. They should probably have used memcpy since they bothered to include extern "C" {} stuff in case it's compiled as C++) \$\endgroup\$ – Peter Cordes Sep 25 at 0:22
  • 3
    \$\begingroup\$ @JL2210 and Rainer:: The other problem with 80-bit long double is that the leading 1 or 0 in the significand (aka mantissa) is explicit, not implied by the exponent field being non-zero. en.wikipedia.org/wiki/… So shifting a bit from the exponent into the mantissa could be a big problem, potentially creating a non-normalized value that isn't a subnormal. I forget what x86/x87 hardware does in that case. If you were using C++, you could use a std::bitset<sizeof(T)*CHAR_BIT> \$\endgroup\$ – Peter Cordes Sep 25 at 0:32
8
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Signed bit-shifting isn't well-defined. I don't think the second half of your answer is valid. \$\endgroup\$ – JL2210 Sep 24 at 17:10
  • \$\begingroup\$ I agree with the top point, but it would be DBL_MANT_DIG-1, not DBL_MANT_DIG. \$\endgroup\$ – JL2210 Sep 24 at 17:33
  • \$\begingroup\$ @JL2210: Is right-shifting a signed bitfield any better defined than a signed int? Doesn't your current code still rely on >> being an arithmetic right shift that shifts in 1s for negative exponents? Yes it's implementation-defined whether it's arithmetic or logical (not undefined). However your code uses GNU C __attribute__((packed)) so you actually can rely on arithmetic right shift if you're intentionally writing in GNU C, not ISO C. gnu.huihoo.org/gcc/gcc-4.9.4/gcc/Integers-implementation.html specifies that "Signed ‘>>’ acts on negative numbers by sign extension." \$\endgroup\$ – Peter Cordes Sep 25 at 1:04
  • 1
    \$\begingroup\$ @JL2210: Yeah, just noticed that. I assumed they would be signed because unsigned right-shift of the 2's complement exponent (after unbiasing) doesn't make sense, does it? A negative exponent (small magnitude float) becomes a very large positive exponent (close-ish to DBL_MAX) after your logical right shift of the exponent field. Or am I missing something and your function does work correctly for inputs like 0.0001? \$\endgroup\$ – Peter Cordes Sep 25 at 1:07
  • 2
    \$\begingroup\$ @JL2210: Hardware that reads the float will effectively undo the bias. It's just an encoding method; the real value is a 2's complement integer. That's what you're working with after undoing the bias (decoding). Yes you have to redo the encoding to make it an IEEE float again, but so what? BTW, if you were to logical right-shift without undoing the bias first, you'd make values twice as close to the smallest possible exponent value instead of twice as close to the middle value (0 after bias). So 1.0f would become 1e-17f or something, if FLT_MIN is about 1e-34f. \$\endgroup\$ – Peter Cordes Sep 25 at 1:22
8
\$\begingroup\$

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.

\$\endgroup\$
  • 1
    \$\begingroup\$ @JL2210: I know you can't have int a:64, that's why I said the underlying type can be wider (than the field) without causing a problem, not vice versa. Apparently my wording is easy to misread, will think about how to fix. \$\endgroup\$ – Peter Cordes Sep 25 at 2:05
  • 1
    \$\begingroup\$ @JL2210: Oh. I meant always fine for the exponent where you were previously using unsigned short. Because both types have the same minimum width in ISO C. I meant that using uint64_t e:8 could give slower-than-necessary code for the exponent. Specifically with MSVC, I was remembering writing this answer: Unresolved external symbol __aullshr when optimization is turned off \$\endgroup\$ – Peter Cordes Sep 25 at 2:09
  • 1
    \$\begingroup\$ @JL2210: Oh yes, I missed that you check separately for == 0.0. I was assuming that you would just use one check on the fast-path to exclude both negative an NaN inputs by requiring that x >= 0.0. That's false for NaN so you can move the isnan check inside an if that won't run for the fast-path. So you can do if (! x>0.0) { further checks to detect +=0.0, NaN, or negative } \$\endgroup\$ – Peter Cordes Sep 25 at 23:28
  • 1
    \$\begingroup\$ @JL2210: Updated my answer to fix that mistake you pointed out, and talk about some of the stuff I've said in comments. And more. \$\endgroup\$ – Peter Cordes Sep 26 at 0:08
  • 1
    \$\begingroup\$ @JL22: Like I said at the top of my answer, sqrt is required by IEEE to be correctly-rounded, just like + - * and /. All x86 sqrt instructions do that. (Except of course rsqrtps which intentionally gives a fast approximation of 1/sqrt(x)). Intel's compiler enables fast-math by default (unlike GCC), but Intel's HW is IEEE-compliant. Perhaps you're thinking of instructions like fsin, which aren't required to be correctly rounded and which in practice have huge relative error near Pi/2. randomascii.wordpress.com/2014/10/09/… \$\endgroup\$ – Peter Cordes Sep 26 at 0:27
4
\$\begingroup\$

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?

\$\endgroup\$
  • 1
    \$\begingroup\$ Under Factoring: That's Newton-Raphson, it detects the amount of error and decreases it. It's not the same as (1-a)/2. \$\endgroup\$ – JL2210 Sep 24 at 15:39
  • 3
    \$\begingroup\$ Soooo. The problem isn't with my answer (and I don't think it should receive your edit, yet); it's with the code you posted. A strict interpretation of the CR policies would suggest that since there's already an answer, you can't edit your code. However, I think instead that you should edit your question to show your actual (non-preprocessed) code, and I'll then amend this answer myself. \$\endgroup\$ – Reinderien Sep 24 at 15:47
  • 2
    \$\begingroup\$ @Reinderien Thanks. I've added a clarification. \$\endgroup\$ – JL2210 Sep 24 at 15:50
  • 1
    \$\begingroup\$ Re. factoring: So my first stab at the math was wrong, and I modified it to reflect reality. \$\endgroup\$ – Reinderien Sep 24 at 16:04
  • 2
    \$\begingroup\$ I immediately know what if(c % 2 == 0) does, but I have to think multiple times to figure out what if (!(c & 1)) does. That's premature optimization and IMHO violates clean code rules (KISS). \$\endgroup\$ – Thomas Weller Sep 26 at 11:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.