flatten 20260225

tjp/core/math/Float16.inl

@@ -0,0 +1,186 @@
+// License: Modified MIT (NON-AI)
+// See the LICENSE file in the root directory for license information.
+// Copyright 2025 Timothy Prepscius
+
+#pragma once
+
+namespace tjp::core {
+
+inline
+Float16::operator float() const
+{
+	return convert(v);
+}
+
+inline
+Float16::operator double() const
+{
+	return (float)*this;
+}
+
+inline
+bool operator <(const Float16 &lhs, const Float16 &rhs)
+{
+	return lhs.ordered() < rhs.ordered();
+}
+
+inline
+bool operator >(const Float16 &lhs, const Float16 &rhs)
+{
+	return lhs.ordered() > rhs.ordered();
+}
+
+inline
+bool operator ==(const Float16 &lhs, const Float16 &rhs)
+{
+	return lhs.v == rhs.v;
+}
+
+inline
+bool operator !=(const Float16 &lhs, const Float16 &rhs)
+{
+	return lhs.v != rhs.v;
+}
+
+// ---
+
+inline
+constexpr Float16 Float16::max()
+{
+	return Float16(Float16::Bits { .v = 0x7BFF });
+}
+
+inline
+constexpr Float16 Float16::negative_max()
+{
+	return Float16(Float16::Bits { .v = 0xFBFF });
+}
+
+inline
+constexpr Float16 Float16::min__()
+{
+	return Float16(Float16::Bits { .v = 0x0400 });
+}
+
+// ---
+
+union Float16_Float32Bits {
+	float f;
+	uint32_t u;
+};
+
+inline
+constexpr float Float16::convert(V v)
+{
+	uint16_t hv = v;
+	uint16_t sign = (hv >> 15) & 0x1;
+	uint16_t exp  = (hv >> 10) & 0x1F;
+	uint16_t frac = hv & 0x3FF;
+
+	uint32_t sign32 = sign << 31;
+	uint32_t exp32 = 0, frac32 = 0;
+
+	if (exp == 0) {
+		if (frac == 0) {
+			// Zero
+			exp32 = 0;
+			frac32 = 0;
+		} else {
+			// Subnormal: normalize
+			int e = -14; // unbiased exponent for half subnormals
+			uint32_t f = frac;
+			while ((f & 0x400) == 0) {
+				f <<= 1;
+				e--;
+			}
+			f &= 0x3FF;
+			exp32 = (e + 127) << 23;
+			frac32 = f << 13;
+		}
+	}
+	else if (exp == 0x1F) {
+		// Inf or NaN
+		exp32 = 0xFF << 23;
+		if (frac == 0) {
+			// Inf
+			frac32 = 0;
+		} else {
+			// NaN – preserve payload, set quiet NaN bit
+			frac32 = (frac << 13) | 0x400000;
+		}
+	}
+	else {
+		// Normalized
+		int e = exp - 15 + 127;
+		exp32 = (uint32_t)e << 23;
+		frac32 = (uint32_t)frac << 13;
+	}
+
+	auto u = sign32 | exp32 | frac32;
+	Float16_Float32Bits result { .u = u };
+	return result.f;
+}
+
+inline
+constexpr Float16::V Float16::convert(float rhs)
+{
+	Float16_Float32Bits result { .f = rhs };
+	uint32_t bits = result.u;
+
+	uint32_t sign = (bits >> 31) & 0x1;
+	int32_t  exp  = (bits >> 23) & 0xFF;
+	uint32_t frac = bits & 0x7FFFFF;
+
+	uint16_t sign16 = sign << 15;
+	uint16_t exp16 = 0;
+	uint16_t frac16 = 0;
+
+	if (exp == 0) {
+		// Zero or subnormal in float
+		exp16 = 0;
+		frac16 = 0; // Flush subnormals to zero
+	}
+	else if (exp == 0xFF) {
+		// Inf or NaN
+		exp16 = 0x1F;
+		frac16 = frac ? 0x200 : 0; // Preserve a quiet NaN bit
+	}
+	else {
+		// Normalized float
+		int newExp = exp - 127 + 15;
+		if (newExp >= 0x1F) {
+			// Overflow → Inf
+			exp16 = 0x1F;
+			frac16 = 0;
+		}
+		else if (newExp <= 0) {
+			// Underflow → subnormal or zero
+			if (newExp < -10) {
+				exp16 = 0;
+				frac16 = 0;
+			}
+			else {
+				// Subnormal half
+				uint32_t mant = (frac | 0x800000) >> (1 - newExp);
+				frac16 = mant >> 13;
+				exp16 = 0;
+			}
+		}
+		else {
+			exp16 = newExp;
+			frac16 = frac >> 13;
+		}
+	}
+
+	Float16::V v = sign16 | (exp16 << 10) | frac16;
+	return v;
+}
+
+inline
+Float16& set(Float16 &lhs, float rhs)
+{
+	lhs.v = Float16::convert(rhs);
+	return lhs;
+}
+
+} // namespace