Teem: air: IEEE 754 floats


teem	/	air	IEEE 754 floats

IEEE 754 is the standard that governs how modern computers represent floating point numbers. One aspect of IEEE 754 which people tend not to exploit is its ability to represent special values like Not-A-Number, and +/-Infinity. The air library has a number of macros and functions which are intended to make working with these quantities easier.

"Existant" values

Parts of teem, especially nrrd, rely on the notion of a floating point value "existing" or not. A value "exists" if it is not a 754 special value, that is, its not infinity, and not NaN. In nrrd, for instance, many operations are defined to only operate on the existant values. Nrrd's average of 1000 numbers, containing one NaN and one infinity, will actually be the average of the 998 "existant" numbers. However, any naive C or C++ implementation will result in the final answer being NaN. Some way of determine if a value exists is required.

Sadly, this requires some care. Ideally, we could define a macro like this:

#define AIR_EXISTS(x) ((x) == (x))

This is supposed to work because x isn't equal to itself if x is a special value, but many compilers will reduce this to true, even at low optimization levels. So, we could try this:

#define AIR_EXISTS(x) (!((x) - (x)))

For special values, x - x is not 0, so this should work. But some compilers (such as 64-bit Irix6 with "-Ofast", and Visual Studio) again decree that AIR_EXISTS() is true for any possible value of x. To be completely thorough, one has to actually look at the bit patterns of the values, at which point you have to be specific to the float/ double type, and correct about their byte-sizes:

#define AIR_EXISTS_F(x) ((*(unsigned int*)&(x) & 0x7f800000) != 0x7f800000)
#define AIR_EXISTS_D(x) (                               \
  (*(airULLong*)&(x) & AIR_ULLONG(0x7ff0000000000000))  \
    != AIR_ULLONG(0x7ff0000000000000))

We want to avoid false negatives for doubles with values greater than FLT_MAX, so if we had to pick one of these to wrap in a single function that was used for both float and double, it would be the second. This is the approach currently taken: on machines for which (!((x) - (x))) doesn't work because of optimizations, we resort to casting the value to double, and looking at its bits (done in airExists_d()). Thus, in air.h we find:

#ifdef WIN32
#define AIR_EXISTS(x) (airExists_d(x))
#else
#define AIR_EXISTS(x) (!((x) - (x)))
#endif

Hopefully, at some soon, teem will be made with the assistance of a highly customized GNU configure script which will determine these things on a per-architecture, per-optimization basis. The above solution works on all platforms that teem compiles on, using the optimization levels that the makefiles supply (generally -O2). A function called airSanity() makes sure that all IEEE 754-related macros (such as AIR_EXISTS()) work correctly.

Not-A-Number

Not-A-Number, or "NaN", doesn't represent a real value, but rather the fact that there is no value to be represented. Various teem libraries, especially nrrd, use NaN as a flag to mean "don't know". This use is completely in keeping with the 754 spec. However, the special semantics of NaN mean that it is not trivial to create a #define or a compile-time constant to hold this value (as is possible with say, M_PI), and it is not trivial to test to see if a value is NaN. This may be why more people don't use NaN as a flag, but air is useful here.

One thing to keep in mind: Section 6.2 of IEEE 754 guarantees that a NaN stored in a float will stay a NaN when stored in a double, and vice versa. Thus, you need not be cautious about the "precision" of your NaN. The same doesn't exactly hold for infinity, though.

How to get NaN: If you want a NaN, so that you can save it into a floating point variable, perhaps as an indicator for "don't know", or "uninitialized", you have a few options.

(fastest) Use AIR_NAN, as in "float nan=AIR_NAN;". AIR_NAN is a #define for airFloatNaN.f, which is the float field of the airFloatNaN union, which is itself a compile-time constant, defined in teem/src/754.c
(slower) Use airFPGen_f(airFP_QNAN), which tells the airFPGen_f() function to generate a float of the class airFP_QNAN. Calling airNaN() is a wrapper around this.
Do it yourself. Store 0.0 in some float variable nan, and then "nan /= nan". You have to create a seperate variable in order to avoid compiler warnings about the division by zero.

How to test for NaN: If you want to see if a variable x of type float is NaN, you can do one of the following. Both of these return non-zero when the passed value is a NaN, and zero when it is not.

(easiest to remember) Use AIR_EXISTS(x), described above. Will be false for NaNs and infinities.
(fastest) Use AIR_ISNAN_F(x), a macro which does bit-masking to look at the exponent and fraction fields of the float
(slower) Use airIsNaN(x). Same as macro, but uses a union to access the exponent and fraction fields.
See if return from int airFPClass_f(float val) is airFP_SNAN or airFP_QNAN.

If you have a variable that is of type double, you must use the second option, because the value will be cast to float as part of the function call airIsNaN, but NOT as part of macro invocation AIR_ISNAN_F.

There are in fact two kinds of NaNs, the signaling NaN and the quiet NaN. Nothing in teem attempts to use the floating point exception handler to work with signaling NaNs, because there isn't a cross-platform C-language API to the floating point hardware. Also, nothing in teem attempts to use the integer value stored in the fraction field of the NaN to represent different kinds of conditions of NaN, so effectively, a NaN is a NaN is a NaN. However, an effort is made to make AIR_NAN be what the local architecture considers a quiet NaN, in case using a signalling NaN could have an impact on performace.

Infinities

A float infinity, when assigned to double, stays infinity (with the same sign), and vice versa. But if a double holds a large value (larger than FLT_MAX) and this is assigned to a float, an infinity will result. Thus, detecting an infinity has to be done either with functions specific to the type of the value, or (being conservative), with functions which take a double argument.

To generate an infinity:

Use AIR_POS_INF or AIR_NEG_INF, a #define for the float component of airFloatPosInf or airFloatNegInf; these are compile-time constant unions.
Use airFPGen_f(airFP_POS_INF) or airFPGen_f(airFP_NEG_INF)

To test if something is infinity:

(easiest to remember) Use AIR_EXISTS(x), described above. Will be false for NaNs and infinities.
Use int airIsInf_f(float f) (or int airIsInf_d(double d)). These return 1 when given a positive infinity, -1 when given a negative infinity, and 0 otherwise.
See if return from int airFPClass_f(float val) (or int airFPClass_d(double val)) is airFP_POS_INF or airFP_NEG_INF.

Determining/generating floating point classes

air.h has an enum for the different classes that a floating point number can belong to:

enum {
  airFP_Unknown,               /*  0: nobody knows */
  airFP_SNAN,                  /*  1: signalling NaN */
  airFP_QNAN,                  /*  2: quiet NaN */
  airFP_POS_INF,               /*  3: positive infinity */
  airFP_NEG_INF,               /*  4: negative infinity */
  airFP_POS_NORM,              /*  5: positive normalized non-zero */
  airFP_NEG_NORM,              /*  6: negative normalized non-zero */
  airFP_POS_DENORM,            /*  7: positive denormalized non-zero */
  airFP_NEG_DENORM,            /*  8: negative denormalized non-zero */
  airFP_POS_ZERO,              /*  9: +0.0, positive zero */
  airFP_NEG_ZERO,              /* 10: -0.0, negative zero */
  airFP_Last                   /* after the last valid one */
};

Values (of type int) from this enum are used in conjunction with:

float airFPGen_f(int cls);
double airFPGen_d(int cls);

These generate float or double values of the given class.

```
int airFPClass_f(float val);
int airFPClass_d(double val);
```
These categorize a given floating-point value into one of the classes.

Ideally, I wouldn't have to reimplement these, since most platforms have them in some form or another. But they are not part of ANSI C, so each platform has its own name for the header file in which these are found, different enum names for the different classes, and even different corresponding numeric values.

Converting between floating-point values and integer triples

There are a number of functions for creating and querying floating point values, written largely for the sake of learning about 754.

```
float airFPPartsToVal_f(int sign, int exp, int frac);
void airFPValToParts_f(int *signP, int *expP, int *fracP, float v);
```
Let you convert between a float and the three integral fields which comprise it: 1-bit sign, 8-bit exponent, and 23-bit fraction.

```
double airFPPartsToVal_d(int sign, int exp, airULLong frac);
void airFPValToParts_d(int *signP, int *expP, airULLong *fracP, double v);
```
Same thing, but for double. airULLong is just another name for unsigned long long, which Microsoft calls an "unsigned __int64". A double is a 1-bit sign, 11-bit exponent, 52-bit fraction.

Printing individual bits

In case you don't trust yourself to read the hex output from the "%x" printf() conversion sequence, you can look at all the individual bits of a float or double:

void airFPFprintf_f(FILE *file, float val);
void airFPFprintf_d(FILE *file, double val);

Sample output from airFPFprintf_f(stdout, 42):

42.000000: class 5; 0x42280000 = sign:0x0, exp:0x84, frac:0x280000 = 
 S < . . Exp . . > < . . . . . . . . . Frac. . . . . . . . . . >
 0 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Sample output from airFPFprintf_d(stdout, 42):

42.000000: class 5; 0x4045000000000000 = 
sign:0x0, exp:0x404, frac:0x5000000000000 = 
S<...Exp...><.......................Frac.......................>
0100000001000101000000000000000000000000000000000000000000000000

Command-line demos are teem/src/air/test/floatprint.c and teem/src/air/test/doubleprint.c