namespace movit {
namespace {
-union fp64 {
- double f;
- unsigned long long ll;
+union fp32 {
+ float f;
+ unsigned int u;
};
-template<class FP16_INT_T,
- int FP16_BIAS, int FP16_MANTISSA_BITS, int FP16_EXPONENT_BITS, int FP16_MAX_EXPONENT,
- int FP64_BIAS, int FP64_MANTISSA_BITS, int FP64_EXPONENT_BITS, int FP64_MAX_EXPONENT>
-inline double fp_upconvert(FP16_INT_T x)
-{
- int sign = x >> (FP16_MANTISSA_BITS + FP16_EXPONENT_BITS);
- int exponent = (x & ((1ULL << (FP16_MANTISSA_BITS + FP16_EXPONENT_BITS)) - 1)) >> FP16_MANTISSA_BITS;
- unsigned long long mantissa = x & ((1ULL << FP16_MANTISSA_BITS) - 1);
-
- int sign64;
- int exponent64;
- unsigned long long mantissa64;
-
- if (exponent == 0) {
- /*
- * Denormals, or zero. Zero is still zero, denormals become
- * ordinary numbers.
- */
- if (mantissa == 0) {
- sign64 = sign;
- exponent64 = 0;
- mantissa64 = 0;
- } else {
- sign64 = sign;
- exponent64 = FP64_BIAS - FP16_BIAS;
- mantissa64 = mantissa << (FP64_MANTISSA_BITS - FP16_MANTISSA_BITS + 1);
-
- /* Normalize the number. */
- while ((mantissa64 & (1ULL << FP64_MANTISSA_BITS)) == 0) {
- --exponent64;
- mantissa64 <<= 1;
- }
-
- /* Clear the now-implicit one-bit. */
- mantissa64 &= ~(1ULL << FP64_MANTISSA_BITS);
- }
- } else if (exponent == FP16_MAX_EXPONENT) {
- /*
- * Infinities or NaN (mantissa=0 => infinity, otherwise NaN).
- * We don't care much about NaNs, so let us just make sure we
- * keep the first bit (which signals signalling/non-signalling
- * in many implementations).
- */
- sign64 = sign;
- exponent64 = FP64_MAX_EXPONENT;
- mantissa64 = mantissa << (FP64_MANTISSA_BITS - FP16_MANTISSA_BITS);
- } else {
- sign64 = sign;
-
- /* Up-conversion is simple. Just re-bias the exponent... */
- exponent64 = exponent + FP64_BIAS - FP16_BIAS;
+} // namespace
- /* ...and convert the mantissa. */
- mantissa64 = mantissa << (FP64_MANTISSA_BITS - FP16_MANTISSA_BITS);
- }
+#ifndef __F16C__
- union fp64 nx;
- nx.ll = ((unsigned long long)sign64 << (FP64_MANTISSA_BITS + FP64_EXPONENT_BITS))
- | ((unsigned long long)exponent64 << FP64_MANTISSA_BITS)
- | mantissa64;
- return nx.f;
-}
-
-unsigned long long shift_right_with_round(unsigned long long x, unsigned shift)
+float fp16_to_fp32(fp16_int_t h)
{
- /* shifts >= 64 need to be special-cased */
- if (shift > 64) {
- return 0;
- } else if (shift == 64) {
- if (x > (1ULL << 63)) {
- return 1;
- } else {
- return 0;
- }
- }
-
- unsigned long long round_part = x & ((1ULL << shift) - 1);
- if (round_part < (1ULL << (shift - 1))) {
- /* round down */
- x >>= shift;
- } else if (round_part > (1ULL << (shift - 1))) {
- /* round up */
- x >>= shift;
- ++x;
+ fp32 magic; magic.u = 113 << 23;
+ unsigned int shifted_exp = 0x7c00 << 13; // exponent mask after shift
+ fp32 o;
+
+ // mantissa+exponent
+ unsigned int shifted = (h.val & 0x7fff) << 13;
+ unsigned int exponent = shifted & shifted_exp;
+
+ // exponent cases
+ o.u = shifted;
+ if (exponent == 0) { // Zero / Denormal
+ o.u += magic.u;
+ o.f -= magic.f;
+ } else if (exponent == shifted_exp) { // Inf/NaN
+ o.u += (255 - 31) << 23;
} else {
- /* round to nearest even number */
- x >>= shift;
- if ((x & 1) != 0) {
- ++x;
- }
+ o.u += (127 - 15) << 23;
}
- return x;
+
+ o.u |= (h.val & 0x8000) << 16; // copy sign bit
+ return o.f;
}
-template<class FP16_INT_T,
- int FP16_BIAS, int FP16_MANTISSA_BITS, int FP16_EXPONENT_BITS, int FP16_MAX_EXPONENT,
- int FP64_BIAS, int FP64_MANTISSA_BITS, int FP64_EXPONENT_BITS, int FP64_MAX_EXPONENT>
-inline FP16_INT_T fp_downconvert(double x)
+fp16_int_t fp32_to_fp16(float x)
{
- union fp64 nx;
- nx.f = x;
- unsigned long long f = nx.ll;
- int sign = f >> (FP64_MANTISSA_BITS + FP64_EXPONENT_BITS);
- int exponent = (f & ((1ULL << (FP64_MANTISSA_BITS + FP64_EXPONENT_BITS)) - 1)) >> FP64_MANTISSA_BITS;
- unsigned long long mantissa = f & ((1ULL << FP64_MANTISSA_BITS) - 1);
-
- int sign16;
- int exponent16;
- unsigned long long mantissa16;
-
- if (exponent == 0) {
- /*
- * Denormals, or zero. The largest possible 64-bit
- * denormal is about +- 2^-1022, and the smallest possible
- * 16-bit denormal is +- 2^-24. Thus, we can safely
- * just set all of these to zero (but keep the sign bit).
- */
- sign16 = sign;
- exponent16 = 0;
- mantissa16 = 0;
- } else if (exponent == FP64_MAX_EXPONENT) {
- /*
- * Infinities or NaN (mantissa=0 => infinity, otherwise NaN).
- * We don't care much about NaNs, so let us just keep the first
- * bit (which signals signalling/ non-signalling) and make sure
- * that we don't coerce NaNs down to infinities.
- */
- if (mantissa == 0) {
- sign16 = sign;
- exponent16 = FP16_MAX_EXPONENT;
- mantissa16 = 0;
- } else {
- sign16 = sign; /* undefined */
- exponent16 = FP16_MAX_EXPONENT;
- mantissa16 = mantissa >> (FP64_MANTISSA_BITS - FP16_MANTISSA_BITS);
- if (mantissa16 == 0) {
- mantissa16 = 1;
- }
- }
- } else {
- /* Re-bias the exponent, and check if we will create a denormal. */
- exponent16 = exponent + FP16_BIAS - FP64_BIAS;
- if (exponent16 <= 0) {
- int shift_amount = FP64_MANTISSA_BITS - FP16_MANTISSA_BITS - exponent16 + 1;
- sign16 = sign;
- exponent16 = 0;
- mantissa16 = shift_right_with_round(mantissa | (1ULL << FP64_MANTISSA_BITS), shift_amount);
-
- /*
- * We could actually have rounded back into the lowest possible non-denormal
- * here, so check for that.
- */
- if (mantissa16 == (1ULL << FP16_MANTISSA_BITS)) {
- exponent16 = 1;
- mantissa16 = 0;
- }
+ fp32 f; f.f = x;
+ unsigned int f32infty = 255 << 23;
+ unsigned int f16max = (127 + 16) << 23;
+ fp32 denorm_magic; denorm_magic.u = ((127 - 15) + (23 - 10) + 1) << 23;
+ unsigned int sign_mask = 0x80000000u;
+ fp16_int_t o = { 0 };
+
+ unsigned int sign = f.u & sign_mask;
+ f.u ^= sign;
+
+ // NOTE all the integer compares in this function can be safely
+ // compiled into signed compares since all operands are below
+ // 0x80000000. Important if you want fast straight SSE2 code
+ // (since there's no unsigned PCMPGTD).
+
+ if (f.u >= f16max) { // result is Inf or NaN (all exponent bits set)
+ o.val = (f.u > f32infty) ? 0x7e00 : 0x7c00; // NaN->qNaN and Inf->Inf
+ } else { // (De)normalized number or zero
+ if (f.u < (113 << 23)) { // resulting FP16 is subnormal or zero
+ // use a magic value to align our 10 mantissa bits at the bottom of
+ // the float. as long as FP addition is round-to-nearest-even this
+ // just works.
+ f.f += denorm_magic.f;
+
+ // and one integer subtract of the bias later, we have our final float!
+ o.val = f.u - denorm_magic.u;
} else {
- /*
- * First, round off the mantissa, since that could change
- * the exponent. We use standard IEEE 754r roundTiesToEven
- * mode.
- */
- sign16 = sign;
- mantissa16 = shift_right_with_round(mantissa, FP64_MANTISSA_BITS - FP16_MANTISSA_BITS);
-
- /* Check if we overflowed and need to increase the exponent. */
- if (mantissa16 == (1ULL << FP16_MANTISSA_BITS)) {
- ++exponent16;
- mantissa16 = 0;
- }
-
- /* Finally, check for overflow, and create +- inf if we need to. */
- if (exponent16 >= FP16_MAX_EXPONENT) {
- exponent16 = FP16_MAX_EXPONENT;
- mantissa16 = 0;
- }
+ unsigned int mant_odd = (f.u >> 13) & 1; // resulting mantissa is odd
+
+ // update exponent, rounding bias part 1
+ f.u += ((15 - 127) << 23) + 0xfff;
+ // rounding bias part 2
+ f.u += mant_odd;
+ // take the bits!
+ o.val = f.u >> 13;
}
}
- return (sign16 << (FP16_MANTISSA_BITS + FP16_EXPONENT_BITS))
- | (exponent16 << FP16_MANTISSA_BITS)
- | mantissa16;
+ o.val |= sign >> 16;
+ return o;
}
-const int FP64_BIAS = 1023;
-const int FP64_MANTISSA_BITS = 52;
-const int FP64_EXPONENT_BITS = 11;
-const int FP64_MAX_EXPONENT = (1 << FP64_EXPONENT_BITS) - 1;
-
-const int FP32_BIAS = 127;
-const int FP32_MANTISSA_BITS = 23;
-const int FP32_EXPONENT_BITS = 8;
-const int FP32_MAX_EXPONENT = (1 << FP32_EXPONENT_BITS) - 1;
-
-const int FP16_BIAS = 15;
-const int FP16_MANTISSA_BITS = 10;
-const int FP16_EXPONENT_BITS = 5;
-const int FP16_MAX_EXPONENT = (1 << FP16_EXPONENT_BITS) - 1;
-
-} // namespace
-
-double fp16_to_fp64(fp16_int_t x)
-{
- return fp_upconvert<fp16_int_t,
- FP16_BIAS, FP16_MANTISSA_BITS, FP16_EXPONENT_BITS, FP16_MAX_EXPONENT,
- FP64_BIAS, FP64_MANTISSA_BITS, FP64_EXPONENT_BITS, FP64_MAX_EXPONENT>(x);
-}
-
-fp16_int_t fp64_to_fp16(double x)
-{
- return fp_downconvert<fp16_int_t,
- FP16_BIAS, FP16_MANTISSA_BITS, FP16_EXPONENT_BITS, FP16_MAX_EXPONENT,
- FP64_BIAS, FP64_MANTISSA_BITS, FP64_EXPONENT_BITS, FP64_MAX_EXPONENT>(x);
-}
-
-double fp32_to_fp64(fp32_int_t x)
-{
- return fp_upconvert<fp32_int_t,
- FP32_BIAS, FP32_MANTISSA_BITS, FP32_EXPONENT_BITS, FP32_MAX_EXPONENT,
- FP64_BIAS, FP64_MANTISSA_BITS, FP64_EXPONENT_BITS, FP64_MAX_EXPONENT>(x);
-}
-
-fp32_int_t fp64_to_fp32(double x)
-{
- return fp_downconvert<fp32_int_t,
- FP32_BIAS, FP32_MANTISSA_BITS, FP32_EXPONENT_BITS, FP32_MAX_EXPONENT,
- FP64_BIAS, FP64_MANTISSA_BITS, FP64_EXPONENT_BITS, FP64_MAX_EXPONENT>(x);
-}
+#endif
} // namespace