git.sesse.net Git - movit/blob - fp16.h

   1 #ifndef _MOVIT_FP16_H
   2 #define _MOVIT_FP16_H 1
   3
   4 #ifdef __F16C__
   5 #include <immintrin.h>
   6 #endif
   7
   8 // Code for converting to and from fp16 (from fp64), without any particular
   9 // machine support, with proper IEEE round-to-even behavior (and correct
  10 // handling of NaNs and infinities). This is needed because some OpenGL
  11 // drivers don't properly round off when asked to convert data themselves.
  12 //
  13 // These routines are originally written by Fabian Giesen, and released by
  14 // him into the public domain;
  15 // see https://fgiesen.wordpress.com/2012/03/28/half-to-float-done-quic/.
  16 // They are quite fast, and can be vectorized if need be; of course, using
  17 // the f16c instructions (see below) will be faster still.
  18
  19 namespace movit {
  20
  21 // structs instead of ints, so that they are not implicitly convertible.
  22 struct fp32_int_t {
  23         unsigned int val;
  24 };
  25 struct fp16_int_t {
  26         unsigned short val;
  27 };
  28
  29 #ifdef __F16C__
  30
  31 // Use the f16c instructions from Haswell if available (and we know that they
  32 // are at compile time).
  33 static inline float fp16_to_fp32(fp16_int_t x)
  34 {
  35         return _cvtsh_ss(x.val);
  36 }
  37
  38 static inline fp16_int_t fp32_to_fp16(float x)
  39 {
  40         fp16_int_t ret;
  41         ret.val = _cvtss_sh(x, 0);
  42         return ret;
  43 }
  44
  45 #else
  46
  47 union fp32 {
  48         float f;
  49         unsigned int u;
  50 };
  51
  52 static inline float fp16_to_fp32(fp16_int_t h)
  53 {
  54         fp32 magic; magic.u = 113 << 23;
  55         unsigned int shifted_exp = 0x7c00 << 13;  // exponent mask after shift
  56         fp32 o;
  57
  58         // mantissa+exponent
  59         unsigned int shifted = (h.val & 0x7fff) << 13;
  60         unsigned int exponent = shifted & shifted_exp;
  61
  62         // exponent cases
  63         o.u = shifted;
  64         if (exponent == 0) {  // Zero / Denormal
  65                 o.u += magic.u;
  66                 o.f -= magic.f;
  67         } else if (exponent == shifted_exp) {  // Inf/NaN
  68                 o.u += (255 - 31) << 23;
  69         } else {
  70                 o.u += (127 - 15) << 23;
  71         }
  72
  73         o.u |= (h.val & 0x8000) << 16;  // copy sign bit
  74         return o.f;
  75 }
  76
  77 static inline fp16_int_t fp32_to_fp16(float x)
  78 {
  79         fp32 f; f.f = x;
  80         unsigned int f32infty = 255 << 23;
  81         unsigned int f16max = (127 + 16) << 23;
  82         fp32 denorm_magic; denorm_magic.u = ((127 - 15) + (23 - 10) + 1) << 23;
  83         unsigned int sign_mask = 0x80000000u;
  84         fp16_int_t o = { 0 };
  85
  86         unsigned int sign = f.u & sign_mask;
  87         f.u ^= sign;
  88
  89         // NOTE all the integer compares in this function can be safely
  90         // compiled into signed compares since all operands are below
  91         // 0x80000000. Important if you want fast straight SSE2 code
  92         // (since there's no unsigned PCMPGTD).
  93
  94         if (f.u >= f16max) {  // result is Inf or NaN (all exponent bits set)
  95                 o.val = (f.u > f32infty) ? 0x7e00 : 0x7c00; // NaN->qNaN and Inf->Inf
  96         } else {  // (De)normalized number or zero
  97                 if (f.u < (113 << 23)) {  // resulting FP16 is subnormal or zero
  98                         // use a magic value to align our 10 mantissa bits at the bottom of
  99                         // the float. as long as FP addition is round-to-nearest-even this
 100                         // just works.
 101                         f.f += denorm_magic.f;
 102
 103                         // and one integer subtract of the bias later, we have our final float!
 104                         o.val = f.u - denorm_magic.u;
 105                 } else {
 106                         unsigned int mant_odd = (f.u >> 13) & 1; // resulting mantissa is odd
 107
 108                         // update exponent, rounding bias part 1
 109                         f.u += (unsigned(15 - 127) << 23) + 0xfff;
 110                         // rounding bias part 2
 111                         f.u += mant_odd;
 112                         // take the bits!
 113                         o.val = f.u >> 13;
 114                 }
 115         }
 116
 117         o.val |= sign >> 16;
 118         return o;
 119 }
 120
 121 #endif
 122
 123 // Overloads for use in templates.
 124 static inline float to_fp32(double x) { return x; }
 125 static inline float to_fp32(float x) { return x; }
 126 static inline float to_fp32(fp16_int_t x) { return fp16_to_fp32(x); }
 127
 128 template<class T> inline T from_fp32(float x);
 129 template<> inline double from_fp32<double>(float x) { return x; }
 130 template<> inline float from_fp32<float>(float x) { return x; }
 131 template<> inline fp16_int_t from_fp32<fp16_int_t>(float x) { return fp32_to_fp16(x); }
 132
 133 template<class From, class To>
 134 inline To convert_float(From x) { return from_fp32<To>(to_fp32(x)); }
 135
 136 template<class Same>
 137 inline Same convert_float(Same x) { return x; }
 138
 139 }  // namespace movit
 140
 141 #endif  // _MOVIT_FP16_H