HalfFloat.h

Go to the documentation of this file.

// Jolt Physics Library (https://github.com/jrouwe/JoltPhysics)
// SPDX-FileCopyrightText: 2021 Jorrit Rouwe
// SPDX-License-Identifier: MIT
 
#pragma once
 
#include <Jolt/Math/Vec4.h>
#include <Jolt/Core/FPException.h>
 
JPH_NAMESPACE_BEGIN
 
using HalfFloat = uint16;
 
// Define half float constant values
static constexpr HalfFloat HALF_FLT_MAX             = 0x7bff;
static constexpr HalfFloat HALF_FLT_MAX_NEGATIVE    = 0xfbff;
static constexpr HalfFloat HALF_FLT_INF             = 0x7c00;
static constexpr HalfFloat HALF_FLT_INF_NEGATIVE    = 0xfc00;
static constexpr HalfFloat HALF_FLT_NANQ            = 0x7e00;
static constexpr HalfFloat HALF_FLT_NANQ_NEGATIVE   = 0xfe00;
 
namespace HalfFloatConversion {
 
// Layout of a float
static constexpr int FLOAT_SIGN_POS = 31;
static constexpr int FLOAT_EXPONENT_POS = 23;
static constexpr int FLOAT_EXPONENT_BITS = 8;
static constexpr int FLOAT_EXPONENT_MASK = (1 << FLOAT_EXPONENT_BITS) - 1;
static constexpr int FLOAT_EXPONENT_BIAS = 127;
static constexpr int FLOAT_MANTISSA_BITS = 23;
static constexpr int FLOAT_MANTISSA_MASK = (1 << FLOAT_MANTISSA_BITS) - 1;
static constexpr int FLOAT_EXPONENT_AND_MANTISSA_MASK = FLOAT_MANTISSA_MASK + (FLOAT_EXPONENT_MASK << FLOAT_EXPONENT_POS);
 
// Layout of half float
static constexpr int HALF_FLT_SIGN_POS = 15;
static constexpr int HALF_FLT_EXPONENT_POS = 10;
static constexpr int HALF_FLT_EXPONENT_BITS = 5;
static constexpr int HALF_FLT_EXPONENT_MASK = (1 << HALF_FLT_EXPONENT_BITS) - 1;
static constexpr int HALF_FLT_EXPONENT_BIAS = 15;
static constexpr int HALF_FLT_MANTISSA_BITS = 10;
static constexpr int HALF_FLT_MANTISSA_MASK = (1 << HALF_FLT_MANTISSA_BITS) - 1;
static constexpr int HALF_FLT_EXPONENT_AND_MANTISSA_MASK = HALF_FLT_MANTISSA_MASK + (HALF_FLT_EXPONENT_MASK << HALF_FLT_EXPONENT_POS);
 
enum ERoundingMode
{
    ROUND_TO_NEG_INF,               
    ROUND_TO_POS_INF,               
    ROUND_TO_NEAREST,               
};
 
template <int RoundingMode>
inline HalfFloat FromFloatFallback(float inV)
{
    // Reinterpret the float as an uint32
    uint32 value = BitCast<uint32>(inV);
 
    // Extract exponent
    uint32 exponent = (value >> FLOAT_EXPONENT_POS) & FLOAT_EXPONENT_MASK;
 
    // Extract mantissa
    uint32 mantissa = value & FLOAT_MANTISSA_MASK;
 
    // Extract the sign and move it into the right spot for the half float (so we can just or it in at the end)
    HalfFloat hf_sign = HalfFloat(value >> (FLOAT_SIGN_POS - HALF_FLT_SIGN_POS)) & (1 << HALF_FLT_SIGN_POS);
 
    // Check NaN or INF
    if (exponent == FLOAT_EXPONENT_MASK) // NaN or INF
        return hf_sign | (mantissa == 0? HALF_FLT_INF : HALF_FLT_NANQ);
 
    // Rebias the exponent for half floats
    int rebiased_exponent = int(exponent) - FLOAT_EXPONENT_BIAS + HALF_FLT_EXPONENT_BIAS;
 
    // Check overflow to infinity
    if (rebiased_exponent >= HALF_FLT_EXPONENT_MASK)
    {
        bool round_up = RoundingMode == ROUND_TO_NEAREST || (hf_sign == 0) == (RoundingMode == ROUND_TO_POS_INF);
        return hf_sign | (round_up? HALF_FLT_INF : HALF_FLT_MAX);
    }
 
    // Check underflow to zero
    if (rebiased_exponent < -HALF_FLT_MANTISSA_BITS)
    {
        bool round_up = RoundingMode != ROUND_TO_NEAREST && (hf_sign == 0) == (RoundingMode == ROUND_TO_POS_INF) && (value & FLOAT_EXPONENT_AND_MANTISSA_MASK) != 0;
        return hf_sign | (round_up? 1 : 0);
    }
 
    HalfFloat hf_exponent;
    int shift;
    if (rebiased_exponent <= 0)
    {
        // Underflow to denormalized number
        hf_exponent = 0;
        mantissa |= 1 << FLOAT_MANTISSA_BITS; // Add the implicit 1 bit to the mantissa
        shift = FLOAT_MANTISSA_BITS - HALF_FLT_MANTISSA_BITS + 1 - rebiased_exponent;
    }
    else
    {
        // Normal half float
        hf_exponent = HalfFloat(rebiased_exponent << HALF_FLT_EXPONENT_POS);
        shift = FLOAT_MANTISSA_BITS - HALF_FLT_MANTISSA_BITS;
    }
 
    // Compose the half float
    HalfFloat hf_mantissa = HalfFloat(mantissa >> shift);
    HalfFloat hf = hf_sign | hf_exponent | hf_mantissa;
 
    // Calculate the remaining bits that we're discarding
    uint remainder = mantissa & ((1 << shift) - 1);
 
    if constexpr (RoundingMode == ROUND_TO_NEAREST)
    {
        // Round to nearest
        uint round_threshold = 1 << (shift - 1);
        if (remainder > round_threshold // Above threshold, we must always round
            || (remainder == round_threshold && (hf_mantissa & 1))) // When equal, round to nearest even
            hf++; // May overflow to infinity
    }
    else
    {
        // Round up or down (truncate) depending on the rounding mode
        bool round_up = (hf_sign == 0) == (RoundingMode == ROUND_TO_POS_INF) && remainder != 0;
        if (round_up)
            hf++; // May overflow to infinity
    }
 
    return hf;
}
 
template <int RoundingMode>
JPH_INLINE HalfFloat FromFloat(float inV)
{
#ifdef JPH_USE_F16C
    FPExceptionDisableOverflow disable_overflow;
    JPH_UNUSED(disable_overflow);
 
    union
    {
        __m128i     u128;
        HalfFloat   u16[8];
    } hf;
    __m128 val = _mm_load_ss(&inV);
    switch (RoundingMode)
    {
    case ROUND_TO_NEG_INF:
        hf.u128 = _mm_cvtps_ph(val, _MM_FROUND_TO_NEG_INF);
        break;
    case ROUND_TO_POS_INF:
        hf.u128 = _mm_cvtps_ph(val, _MM_FROUND_TO_POS_INF);
        break;
    case ROUND_TO_NEAREST:
        hf.u128 = _mm_cvtps_ph(val, _MM_FROUND_TO_NEAREST_INT);
        break;
    }
    return hf.u16[0];
#else
    return FromFloatFallback<RoundingMode>(inV);
#endif
}
 
inline Vec4 ToFloatFallback(UVec4Arg inValue)
{
    // Unpack half floats to 4 uint32's
    UVec4 value = inValue.Expand4Uint16Lo();
 
    // Normal half float path, extract the exponent and mantissa, shift them into place and update the exponent bias
    UVec4 exponent_mantissa = UVec4::sAnd(value, UVec4::sReplicate(HALF_FLT_EXPONENT_AND_MANTISSA_MASK)).LogicalShiftLeft<FLOAT_EXPONENT_POS - HALF_FLT_EXPONENT_POS>() + UVec4::sReplicate((FLOAT_EXPONENT_BIAS - HALF_FLT_EXPONENT_BIAS) << FLOAT_EXPONENT_POS);
 
    // Denormalized half float path, renormalize the float
    UVec4 exponent_mantissa_denormalized = ((exponent_mantissa + UVec4::sReplicate(1 << FLOAT_EXPONENT_POS)).ReinterpretAsFloat() - UVec4::sReplicate((FLOAT_EXPONENT_BIAS - HALF_FLT_EXPONENT_BIAS + 1) << FLOAT_EXPONENT_POS).ReinterpretAsFloat()).ReinterpretAsInt();
 
    // NaN / INF path, set all exponent bits
    UVec4 exponent_mantissa_nan_inf = UVec4::sOr(exponent_mantissa, UVec4::sReplicate(FLOAT_EXPONENT_MASK << FLOAT_EXPONENT_POS));
 
    // Get the exponent to determine which of the paths we should take
    UVec4 exponent_mask = UVec4::sReplicate(HALF_FLT_EXPONENT_MASK << HALF_FLT_EXPONENT_POS);
    UVec4 exponent = UVec4::sAnd(value, exponent_mask);
    UVec4 is_denormalized = UVec4::sEquals(exponent, UVec4::sZero());
    UVec4 is_nan_inf = UVec4::sEquals(exponent, exponent_mask);
 
    // Select the correct result
    UVec4 result_exponent_mantissa = UVec4::sSelect(UVec4::sSelect(exponent_mantissa, exponent_mantissa_nan_inf, is_nan_inf), exponent_mantissa_denormalized, is_denormalized);
 
    // Extract the sign bit and shift it to the left
    UVec4 sign = UVec4::sAnd(value, UVec4::sReplicate(1 << HALF_FLT_SIGN_POS)).LogicalShiftLeft<FLOAT_SIGN_POS - HALF_FLT_SIGN_POS>();
 
    // Construct the float
    return UVec4::sOr(sign, result_exponent_mantissa).ReinterpretAsFloat();
}
 
JPH_INLINE Vec4 ToFloat(UVec4Arg inValue)
{
#if defined(JPH_USE_F16C)
    return _mm_cvtph_ps(inValue.mValue);
#elif defined(JPH_USE_NEON)
    return vcvt_f32_f16(vreinterpret_f16_u32(vget_low_u32(inValue.mValue)));
#else
    return ToFloatFallback(inValue);
#endif
}
 
} // HalfFloatConversion
 
JPH_NAMESPACE_END