Move some fp16 math utils here.

math/math_util.h

@@ -114,7 +114,72 @@ inline float round_ieee_754(float d) {
 	return i + 1.0f;
 }
 
+union FP32 {
+	uint32_t u;
+	float f;
+};
 
+struct FP16 {
+	uint16_t u;
+};
+
+// magic code from ryg: http://fgiesen.wordpress.com/2012/03/28/half-to-float-done-quic/
+// See also SSE2 version: https://gist.github.com/rygorous/2144712
+inline FP32 half_to_float_fast5(FP16 h)
+{
+	static const FP32 magic = { (127 + (127 - 15)) << 23 };
+	static const FP32 was_infnan = { (127 + 16) << 23 };
+	FP32 o;
+	o.u = (h.u & 0x7fff) << 13;     // exponent/mantissa bits
+	o.f *= magic.f;                 // exponent adjust
+	if (o.f >= was_infnan.f)        // make sure Inf/NaN survive
+		o.u |= 255 << 23;
+	o.u |= (h.u & 0x8000) << 16;    // sign bit
+	return o;
+}
+
+inline float ExpandHalf(uint16_t half) {
+	FP16 fp16;
+	fp16.u = half;
+	FP32 fp = half_to_float_fast5(fp16);
+	return fp.f;
+}
+
+// More magic code: https://gist.github.com/rygorous/2156668
+inline FP16 float_to_half_fast3(FP32 f)
+{
+	static const FP32 f32infty = { 255 << 23 };
+	static const FP32 f16infty = { 31 << 23 };
+	static const FP32 magic = { 15 << 23 };
+	static const uint32_t sign_mask = 0x80000000u;
+	static const uint32_t round_mask = ~0xfffu;
+	FP16 o = { 0 };
+
+	uint32_t sign = f.u & sign_mask;
+	f.u ^= sign;
+
+	if (f.u >= f32infty.u) // Inf or NaN (all exponent bits set)
+		o.u = (f.u > f32infty.u) ? 0x7e00 : 0x7c00; // NaN->qNaN and Inf->Inf
+	else // (De)normalized number or zero
+	{
+		f.u &= round_mask;
+		f.f *= magic.f;
+		f.u -= round_mask;
+		if (f.u > f16infty.u) f.u = f16infty.u; // Clamp to signed infinity if overflowed
+
+		o.u = f.u >> 13; // Take the bits!
+	}
+
+	o.u |= sign >> 16;
+	return o;
+}
+
+inline uint16_t ShrinkToHalf(float full) {
+	FP32 fp32;
+	fp32.f = full;
+	FP16 fp = float_to_half_fast3(fp32);
+	return fp.u;
+}
 
 // FPU control.
 void EnableFZ();