Appendix B Converting float to half

A Case Study in Rounding

Studying the float→half conversion operation is a useful way to learn the details of floating point encodings and rounding. Because it’s a simple unary operation, we can focus on the encoding and rounding without getting distracted by the details of floating point arithmetic and the precision of intermediate representations.

When converting from float to half, the correct output for any float too large to represent is half infinity; any float too small to represent as a half (even a denormal half) must be clamped to 0.0. The maximum float that rounds to half 0.0 is 0x32FFFFFF, or 2.98×10^-8, while the smallest float that rounds to half infinity is 65520.0. float values inside this range can be converted to half by propagating the sign bit, rebiasing the exponent (since float has an 8-bit exponent biased by 127 and half has a 5-bit exponent biased by 15) and rounding the float mantissa to the nearest half mantissa value. Rounding is straightforward in all cases except when the input value falls exactly between the two possible output values; when this case applies, the IEEE standard specifies rounding to the “nearest even” value. In decimal arithmetic, this would mean rounding 1.5 to 2.0, but also rounding 2.5 to 2.0 and (for example) rounding 0.5 to 0.0.

Listing B-1 shows a C routine that exactly replicates the float-to-half conversion operation. It is built around a set of symbolic constants that name the exponent shifts, bit counts, biases, and masks of the two formats, so the shifts and comparisons in the body read in terms of the formats’ fields rather than hard-coded magic numbers. The routine is adapted from a formulation Norbert Juffa posted to the NVIDIA Developer Forums, with the hard-coded constants replaced by the symbolic ones used here; Juffa’s original is reproduced under its BSD 2-Clause license in the note that follows.¹

The macro LG_MAKE_MASK, used in Listing B-1, creates a mask with a given bit count:

#define LG_MAKE_MASK(bits) ((1<<(bits))-1)

A volatile union is used to treat the same 32-bit value as float and unsigned int; idioms such as *((float *) (&u)) are not portable.

The routine copies the input’s sign bit into the output, then dispatches on the input’s exponent field. If every exponent bit is set, the input is INF or NaN: infinity is propagated (preserving the sign already placed in the output), and anything else becomes a canonical NaN. Otherwise, if the input’s biased exponent falls below f32MinRNonzero – the smallest float exponent whose values can round to a nonzero half, halfway between 0.0 and the smallest half denormal – the input is too small to represent, and the result is the signed 0.0 already sitting in the output. This single threshold replaces the pair of clamps that a more literal implementation would apply at both ends of the representable range.

For the values that remain, the routine makes the input mantissa’s implicit leading 1 explicit and computes shift, the unbiased input exponent. If shift exceeds f16MaxExp, the value overflows to half infinity. Otherwise, the normal and denormal cases must be handled in distinct ways. For a normal result, the mantissa is shifted right by f32ExpShift − f16ExpShift and the rebiased exponent (f16ExpBias − 1 + shift) is added into the exponent field; the bias is reduced by one because the explicit leading 1, shifted up into the exponent field, supplies the missing count. For a denormal result (shift < f16MinExp), the mantissa is shifted right by an additional f16MinExp − shift places, so the implicit 1 lands at the correct denormal position.

In both cases, the bits shifted off the bottom of the mantissa are retained, MSB-aligned, in the working value. Rounding then reduces to a single comparison: 0x80000000 is exactly half an ULP of the result, so if the residual exceeds it, the result is rounded up; if it equals that value exactly – a true tie – the result is rounded up only when the output’s least significant bit is already set. This is the IEEE round-to-nearest-even rule, and expressing it as a comparison against the half-ULP point sidesteps the separate rounding-mask construction that earlier versions of this routine used.

// Make mask out of bit count#define LG_MAKE_MASK(bits) ((1<<(bits))-1) const int f16ExpShift = 10;const int f16MantissaBits = 10;const int f16ExpBits = 5; const int f16ExpBias = 15;const int f16MinExp = -14;const int f16MaxExp = 15;const int f16SignMask = 0x8000; const int f32ExpShift = 23;const int f32MantissaBits = 23;const int f32ExpBits = 8;const int f32ExpBias = 127;const int f32SignMask = 0x80000000; // All bits of the exponent field set (for float16, also infinity).const int f16ExpMask = LG_MAKE_MASK(f16ExpBits) << f16ExpShift;const int f32ExpMask = LG_MAKE_MASK(f32ExpBits) << f32ExpShift; // Smallest float32 exponent field whose values may round to a// nonzero float16 (halfway between 0 and the smallest denormal)const int f32MinRNonzero = (f16MinExp-f16MantissaBits-1+f32ExpBias) <<    f32ExpShift; unsigned shortConvertFloatToHalf( float f ){    // Use a volatile union to portably coerce the 32-bit    // float into a 32-bit integer.    volatile union {        float f;        unsigned int u;    } uf;    uf.f = f;     unsigned int ia = uf.u;    unsigned short ir;     // start by propagating the sign bit    ir = (ia >> 16) & f16SignMask;     if ( (ia & f32ExpMask) == f32ExpMask ) {    // INF or NaN        if ( (ia & ~f32SignMask) == f32ExpMask ) {            ir |= f16ExpMask;                   // INF - propagate sign        }        else {            ir = f16ExpMask | LG_MAKE_MASK(f16MantissaBits);  // canonical NaN        }    }    else if ( (ia & f32ExpMask) >= (unsigned int) f32MinRNonzero ) {        int shift = (int) ((ia >> f32ExpShift) &            LG_MAKE_MASK(f32ExpBits)) - f32ExpBias;        if ( shift > f16MaxExp ) {            ir |= f16ExpMask;                   // overflow - round to infinity        }        else {            // extract mantissa and make the implicit 1 explicit            ia = (ia & LG_MAKE_MASK(f32MantissaBits)) | (1<<f32ExpShift);            if ( shift < f16MinExp ) {          // Case 1: float16 denormal                int RelativeShift = f32ExpShift-f16ExpShift+f16MinExp-shift;                ir |= ia >> RelativeShift;                ia = ia << (32 - RelativeShift);            }            else {                              // Case 2: float16 normal                int RelativeShift = f32ExpShift-f16ExpShift;                ir |= ia >> RelativeShift;                ia = ia << (32 - RelativeShift);                // f16ExpBias-1, since the explicit 1 shifted into                // the exponent field adds 1 to the exponent                ir = ir + ((f16ExpBias-1+shift) << f16ExpShift);            }            // Round to nearest even: ia holds the shifted-out bits,            // MSB-aligned; 0x80000000 is exactly half an ULP of the result.            if ( (ia > 0x80000000) || ((ia == 0x80000000) && (ir & 1)) ) {                ir++;            }        }    }    return ir;}

In practice, developers should convert float to half by using the __float2half() intrinsic. This sample routine is provided purely to aid in understanding floating point layout and rounding; also, examining all the special-case code for INF, NaN and denormal values helps to illustrate why these features of the IEEE spec have been controversial since its inception: they make hardware slower, more costly, or both due to increased silicon area and engineering effort for validation.

In the code accompanying this book, the ConvertFloatToHalf() routine of Listing B-1 is incorporated into a program called float_to_float16.cu that tests its output for every 32-bit floating point value.