GF16: the primary format

You will learn

The 16-bit GoldenFloat the family marks primary: 1 + 6 + 9, bias 31, and its special codes.

Try it

Find SIGN_MASK, EXP_MASK and MANT_MASK in gf16.t27 and check that together they cover all 16 bits. Then find the code of NaN and say which field makes it not an infinity.

specs/numeric/gf16.t27

// SPDX-License-Identifier: Apache-2.0
; gf16.t27 — GoldenFloat16 Encode/Decode
; GF16: 16-bit floating point with 1 sign + 6 exponent + 9 mantissa
; Bit layout: [S(1) E(6) M(9)] = [15:15][14:9][8:0]
; φ² + 1/φ² = 3 | TRINITY

module triformat-gf16;

// ============================================================================
// Constants
// ============================================================================

pub const SIGN_SHIFT : u8 = 15;
pub const EXP_SHIFT : u8 = 9;
pub const MANT_SHIFT : u8 = 0;

pub const SIGN_MASK : u16 = 0x8000;   // 1 << 15
pub const EXP_MASK : u16 = 0x7E00;    // 0b111111 << 9
pub const MANT_MASK : u16 = 0x01FF;   // 0b111111111

pub const EXP_MAX : u8 = 0x3F;       // 63 (all ones in 6 bits)
pub const EXP_MIN : u8 = 0x00;

pub const BIAS : i8 = 31;             // Exponent bias for GF16
pub const SPECIAL_EXP : u8 = 0x3F;    // All ones = special (Inf/NaN)

pub const MANT_DIVISOR : u16 = 512;  // 2^9
pub const MANT_DIVISOR_SHIFT : u8 = 9; // log2(512)

pub const PHI_BIAS : u16 = 60;        // Phi-optimized rounding bias

// GF16 special values
pub const GF16_ZERO_POS : u16 = 0x0000;
pub const GF16_ZERO_NEG : u16 = 0x8000;
pub const GF16_INF_POS : u16 = 0x7E00;
pub const GF16_INF_NEG : u16 = 0xFE00;
pub const GF16_NAN : u16 = 0xFE01;    // Sign + all exp + mantissa != 0

// ============================================================================
// Types
// ============================================================================

pub const GF16 = u16;

// ============================================================================
// Lookup Tables
// ============================================================================

// Powers of 2 for exponents 0-31
pub const pow2_table : [32]u16 = [32]u16{
    0x3C00, 0x3D00, 0x3D80, 0x3E00, 0x3E40, 0x3E80, 0x3EC0, 0x3F00,
    0x3F40, 0x3F80, 0x3FC0, 0x3FE0, 0x3FF0, 0x4000, 0x4040, 0x4080,
    0x40C0, 0x4100, 0x4140, 0x4180, 0x41C0, 0x4200, 0x4240, 0x4280,
    0x42C0, 0x4300, 0x4340, 0x4380, 0x43C0, 0x4400, 0x4440, 0x4480,
};

// ============================================================================
// Functions
// ============================================================================

// gf16_extract_sign(gf16: GF16) → i8
// Extract sign bit (bit 15)
// Returns: 0 for positive, -1 for negative
pub fn gf16_extract_sign(gf16: GF16) i8 {
    const bit = (gf16 >> SIGN_SHIFT) & 1;
    return if (bit != 0) -1 else 0;
}

// gf16_extract_exponent(gf16: GF16) → i8
// Extract exponent bits (bits 14-9)
// Returns: 0-63
pub fn gf16_extract_exponent(gf16: GF16) i8 {
    return @as(i8, @intCast((gf16 >> EXP_SHIFT) & EXP_MASK));
}

// gf16_extract_mantissa(gf16: GF16) → i16
// Extract mantissa bits (bits 8-0)
// Returns: 0-511
pub fn gf16_extract_mantissa(gf16: GF16) i16 {
    return @as(i16, gf16 & MANT_MASK);
}

// gf16_from_components(sign: i8, exp: i8, mant: i16) → GF16
// Assemble GF16 from sign, exponent, mantissa
pub fn gf16_from_components(sign: i8, exp: i8, mant: i16) GF16 {
    const sign_bit = if (sign < 0) 1 else 0;
    return (@as(GF16, @intCast(sign_bit)) << SIGN_SHIFT) |
           (@as(GF16, @intCast(exp)) << EXP_SHIFT) |
           @as(GF16, @intCast(mant));
}

// gf16_is_zero(gf16: GF16) → bool
// Check if GF16 is zero (positive or negative)
pub fn gf16_is_zero(gf16: GF16) bool {
    return gf16 == GF16_ZERO_POS or gf16 == GF16_ZERO_NEG;
}

// gf16_is_special(gf16: GF16) → bool
// Check if GF16 is Inf or NaN (exp == 63)
pub fn gf16_is_special(gf16: GF16) bool {
    return gf16_extract_exponent(gf16) == EXP_MAX;
}

// gf16_encode_f32(f32: f32) → GF16
// Encode IEEE 754 single precision to GF16
// Round-to-nearest, ties to even
// Range: 2^-31 to 2^32 (normal), subnormals flushed to zero
pub fn gf16_encode_f32(value: f32) GF16 {
    // Handle zero
    if (value == 0.0) {
        return if (std.math.signbit(value)) GF16_ZERO_NEG else GF16_ZERO_POS;
    }

    // Extract sign
    const sign = if (value < 0.0) -1 else 0;
    const abs_value = if (value < 0.0) -value else value;

    // Get f32 components
    const f32_bits: u32 = @bitCast(abs_value);
    var f32_exp: i8 = @intCast((f32_bits >> 23) & 0xFF) - 127;
    var f32_mant: u32 = f32_bits & 0x7FFFFF;

    // Convert exp from f32 bias (127) to GF16 bias (31)
    // gf16_exp = f32_exp + 31 - 127 = f32_exp - 96
    var gf16_exp: u8 = @as(u8, @intCast(f32_exp - 96));

    // Clamp exponent
    if (gf16_exp < 0) {
        gf16_exp = 0;  // Underflow to zero
    } else if (gf16_exp > EXP_MAX) {
        gf16_exp = EXP_MAX;  // Overflow to Inf
    }

    // Extract mantissa and scale to 9 bits
    // f32 mantissa is 23 bits, GF16 needs 9 bits
    // Shift right by 14 bits (23 - 9 = 14)
    var mant = @as(u16, @intCast(f32_mant >> 14));

    // Round-to-nearest
    const discarded = f32_mant & 0x3FFF;
    if ((discarded & 0x2000) != 0) {
        mant += 1;
        if (mant > MANT_MASK) {
            mant = 0;
            if (gf16_exp < EXP_MAX) {
                gf16_exp += 1;
            }
        }
    }

    return gf16_from_components(sign, gf16_exp, mant);
}

// gf16_decode_to_f32(gf16: GF16) → f32
// Decode GF16 to IEEE 754 single precision
pub fn gf16_decode_to_f32(gf16: GF16) f32 {
    // Handle zero
    if (gf16_is_zero(gf16)) {
        const sign = gf16_extract_sign(gf16);
        return if (sign < 0) -0.0 else 0.0;
    }

    // Handle special values (Inf/NaN)
    if (gf16_is_special(gf16)) {
        const mant = gf16_extract_mantissa(gf16);
        const sign = gf16_extract_sign(gf16);
        if (mant == 0) {
            // Infinity
            return if (sign < 0) -std.math.inf(f32) else std.math.inf(f32);
        } else {
            // NaN
            return std.math.nan(f32);
        }
    }

    // Normal number: value = (-1)^s * (1 + m/2^9) * 2^(e - 31)
    const sign = gf16_extract_sign(gf16);
    const exp = gf16_extract_exponent(gf16);
    const mant = gf16_extract_mantissa(gf16);

    const sign_mult = if (sign < 0) -1.0 else 1.0;
    const mant_mult = 1.0 + @as(f32, @floatFromInt(mant)) / 512.0;
    const exp_mult = @as(f32, @exp2(f32, @floatFromInt(exp - BIAS)));

    return sign_mult * mant_mult * exp_mult;
}

// gf16_round_phi(value: f32) → GF16
// Phi-optimized rounding for GF16
// Uses golden ratio bias for rounding decisions instead of standard round-to-nearest
// Bias = (1/φ - 0.5) * scale, where 1/φ ≈ 0.618
// This improves numerical stability for sacred physics calculations
pub fn gf16_round_phi(value: f32) GF16 {
    // Handle zero
    if (value == 0.0) {
        return if (std.math.signbit(value)) GF16_ZERO_NEG else GF16_ZERO_POS;
    }

    // Extract sign
    const sign = if (value < 0.0) -1 else 0;
    const abs_value = if (value < 0.0) -value else value;

    // Get f32 components
    const f32_bits: u32 = @bitCast(abs_value);
    var f32_exp: i8 = @intCast((f32_bits >> 23) & 0xFF) - 127;
    var f32_mant: u32 = f32_bits & 0x7FFFFF;

    // Convert exp from f32 bias (127) to GF16 bias (31)
    var gf16_exp: u8 = @as(u8, @intCast(f32_exp - 96));

    // Clamp exponent
    if (gf16_exp < 0) {
        gf16_exp = 0;
    } else if (gf16_exp > EXP_MAX) {
        gf16_exp = EXP_MAX;
    }

    // Add implied 1 for normalization
    const normalized_mant: u32 = f32_mant | 0x00800000;

    // Scale to 9 bits with phi bias
    var mant = @as(u16, @intCast((normalized_mant >> 15) + PHI_BIAS));

    // Check for overflow and adjust
    if (mant > MANT_MASK) {
        mant = 0;
        if (gf16_exp < EXP_MAX) {
            gf16_exp += 1;
        } else {
            gf16_exp = EXP_MAX;  // Overflow to Inf
        }
    }

    return gf16_from_components(sign, gf16_exp, mant);
}

// gf16_is_inf(gf16: GF16) → bool
// Check if GF16 represents infinity
pub fn gf16_is_inf(gf16: GF16) bool {
    const exp = gf16_extract_exponent(gf16);
    const mant = gf16_extract_mantissa(gf16);
    return (exp == EXP_MAX) and (mant == 0);
}

// gf16_is_nan(gf16: GF16) → bool
// Check if GF16 represents NaN (Not a Number)
pub fn gf16_is_nan(gf16: GF16) bool {
    const exp = gf16_extract_exponent(gf16);
    const mant = gf16_extract_mantissa(gf16);
    return (exp == EXP_MAX) and (mant != 0);
}

// gf16_is_negative(gf16: GF16) → bool
// Check if GF16 is negative (excluding negative zero)
pub fn gf16_is_negative(gf16: GF16) bool {
    const sign = gf16_extract_sign(gf16);
    return (sign < 0) and !gf16_is_zero(gf16);
}

// gf16_is_positive(gf16: GF16) → bool
// Check if GF16 is positive (excluding positive zero)
pub fn gf16_is_positive(gf16: GF16) bool {
    const sign = gf16_extract_sign(gf16);
    return (sign >= 0) and !gf16_is_zero(gf16);
}

// gf16_negate(gf16: GF16) → GF16
// Negate a GF16 value (flip sign bit)
pub fn gf16_negate(gf16: GF16) GF16 {
    return gf16 ^ SIGN_MASK;
}

// gf16_abs(gf16: GF16) → GF16
// Absolute value of GF16 (clear sign bit)
pub fn gf16_abs(gf16: GF16) GF16 {
    return gf16 & ~SIGN_MASK;
}

// gf16_copy_sign(gf16: GF16, sign_source: GF16) → GF16
// Copy sign from sign_source to gf16 value
pub fn gf16_copy_sign(gf16: GF16, sign_source: GF16) GF16 {
    const sign_mask = sign_source & SIGN_MASK;
    const value_mask = gf16 & ~SIGN_MASK;
    return value_mask | sign_mask;
}

// gf16_max(a: GF16, b: GF16) → GF16
// Return the greater of two GF16 values
pub fn gf16_max(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a)) return b;
    if (gf16_is_nan(b)) return a;

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);

    if (a_val >= b_val) return a else return b;
}

// gf16_min(a: GF16, b: GF16) → GF16
// Return the smaller of two GF16 values
pub fn gf16_min(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a)) return b;
    if (gf16_is_nan(b)) return a;

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);

    if (a_val <= b_val) return a else return b;
}

// gf16_add(a: GF16, b: GF16) → GF16
// Add two GF16 values (decode, add, re-encode)
// Returns NaN if either operand is NaN, Inf if overflow
pub fn gf16_add(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return GF16_NAN;
    if (gf16_is_inf(a) and gf16_is_inf(b)) {
        // Inf + Inf = NaN (if same sign)
        // Inf + (-Inf) = NaN
        const a_sign = gf16_extract_sign(a);
        const b_sign = gf16_extract_sign(b);
        return if (a_sign == b_sign) a else GF16_NAN;
    }
    if (gf16_is_inf(a)) return a;
    if (gf16_is_inf(b)) return b;

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const result = a_val + b_val;

    return gf16_encode_f32(result);
}

// gf16_sub(a: GF16, b: GF16) → GF16
// Subtract two GF16 values (decode, subtract, re-encode)
// Returns NaN if either operand is NaN, Inf if overflow
pub fn gf16_sub(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return GF16_NAN;
    if (gf16_is_inf(a) and gf16_is_inf(b)) {
        // Inf - Inf = NaN
        return GF16_NAN;
    }
    if (gf16_is_inf(a)) return a;
    if (gf16_is_inf(b)) {
        // -Inf + something = Inf (with sign flip)
        return gf16_negate(b);
    }

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const result = a_val - b_val;

    return gf16_encode_f32(result);
}

// gf16_mul(a: GF16, b: GF16) → GF16
// Multiply two GF16 values (decode, multiply, re-encode)
// Returns NaN if either operand is NaN
pub fn gf16_mul(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return GF16_NAN;
    if (gf16_is_zero(a) or gf16_is_zero(b)) {
        // 0 * x = 0, with sign handling
        const a_sign = gf16_extract_sign(a);
        const b_sign = gf16_extract_sign(b);
        const result_sign = a_sign ^ b_sign;
        return if (result_sign != 0) GF16_ZERO_NEG else GF16_ZERO_POS;
    }
    if (gf16_is_inf(a) or gf16_is_inf(b)) {
        // Inf * non-zero = Inf, with sign handling
        const a_sign = gf16_extract_sign(a);
        const b_sign = gf16_extract_sign(b);
        const result_sign = a_sign ^ b_sign;
        return if (result_sign != 0) GF16_INF_NEG else GF16_INF_POS;
    }

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const result = a_val * b_val;

    return gf16_encode_f32(result);
}

// gf16_div(a: GF16, b: GF16) → GF16
// Divide two GF16 values (decode, divide, re-encode)
// Returns NaN if division by zero or either operand is NaN
// Returns Inf if numerator is Inf and denominator is finite non-zero
pub fn gf16_div(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return GF16_NAN;
    if (gf16_is_zero(b)) {
        // Division by zero = Inf with sign of a
        const a_sign = gf16_extract_sign(a);
        return if (a_sign != 0) GF16_INF_NEG else GF16_INF_POS;
    }
    if (gf16_is_inf(a)) {
        // Inf / finite = Inf, with sign handling
        const a_sign = gf16_extract_sign(a);
        const b_sign = gf16_extract_sign(b);
        const result_sign = a_sign ^ b_sign;
        return if (result_sign != 0) GF16_INF_NEG else GF16_INF_POS;
    }
    if (gf16_is_inf(b)) {
        // Finite / Inf = 0, with sign handling
        const a_sign = gf16_extract_sign(a);
        const b_sign = gf16_extract_sign(b);
        const result_sign = a_sign ^ b_sign;
        return if (result_sign != 0) GF16_ZERO_NEG else GF16_ZERO_POS;
    }

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const result = a_val / b_val;

    return gf16_encode_f32(result);
}

// gf16_fma(a: GF16, b: GF16, c: GF16) → GF16
// Fused multiply-add: a * b + c with single rounding
// More accurate than separate mul and add
pub fn gf16_fma(a: GF16, b: GF16, c: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b) or gf16_is_nan(c)) return GF16_NAN;

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const c_val = gf16_decode_to_f32(c);

    // Handle special cases
    if (gf16_is_zero(a) or gf16_is_zero(b)) {
        return gf16_add(c, gf16_encode_f32(0.0));
    }

    // Compute a * b + c
    const product = a_val * b_val;
    const result = product + c_val;

    return gf16_encode_f32(result);
}

// gf16_sqrt(a: GF16) → GF16
// Square root of GF16 value
// Returns NaN for negative values, Inf for infinity
pub fn gf16_sqrt(a: GF16) GF16 {
    if (gf16_is_nan(a)) return GF16_NAN;
    if (gf16_is_inf(a) and !gf16_is_negative(a)) return a;
    if (gf16_is_inf(a)) return GF16_NAN; // -Inf sqrt = NaN
    if (gf16_is_zero(a)) return a;
    if (gf16_is_negative(a)) return GF16_NAN;

    const a_val = gf16_decode_to_f32(a);
    const result = @sqrt(a_val);

    return gf16_encode_f32(result);
}

// gf16_square(a: GF16) → GF16
// Square of GF16 value
// Uses gf16_mul internally
pub fn gf16_square(a: GF16) GF16 {
    return gf16_mul(a, a);
}

// gf16_eq(a: GF16, b: GF16) → bool
// Equality comparison for GF16
// NaN values are never equal to anything (including themselves)
pub fn gf16_eq(a: GF16, b: GF16) bool {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return false;
    // For zero values, treat +0 and -0 as equal
    if (gf16_is_zero(a) and gf16_is_zero(b)) return true;
    return a == b;
}

// gf16_ne(a: GF16, b: GF16) → bool
// Not-equal comparison for GF16
// NaN values are not equal to anything (including themselves)
pub fn gf16_ne(a: GF16, b: GF16) bool {
    return !gf16_eq(a, b);
}

// gf16_lt(a: GF16, b: GF16) → bool
// Less-than comparison for GF16
// Returns false if either operand is NaN
pub fn gf16_lt(a: GF16, b: GF16) bool {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return false;
    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    return a_val < b_val;
}

// gf16_le(a: GF16, b: GF16) → bool
// Less-than-or-equal comparison for GF16
// Returns false if either operand is NaN
pub fn gf16_le(a: GF16, b: GF16) bool {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return false;
    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    return a_val <= b_val;
}

// gf16_gt(a: GF16, b: GF16) → bool
// Greater-than comparison for GF16
// Returns false if either operand is NaN
pub fn gf16_gt(a: GF16, b: GF16) bool {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return false;
    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    return a_val > b_val;
}

// gf16_ge(a: GF16, b: GF16) → bool
// Greater-than-or-equal comparison for GF16
// Returns false if either operand is NaN
pub fn gf16_ge(a: GF16, b: GF16) bool {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return false;
    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    return a_val >= b_val;
}

// gf16_floor(a: GF16) → GF16
// Round down to the nearest integer (toward -inf)
// Returns NaN for NaN input, unchanged for Inf/-Inf
pub fn gf16_floor(a: GF16) GF16 {
    if (gf16_is_nan(a)) return GF16_NAN;
    if (gf16_is_inf(a)) return a;
    if (gf16_is_zero(a)) return a;

    const a_val = gf16_decode_to_f32(a);
    const result = @floor(a_val);

    return gf16_encode_f32(result);
}

// gf16_ceil(a: GF16) → GF16
// Round up to the nearest integer (toward +inf)
// Returns NaN for NaN input, unchanged for Inf/-Inf
pub fn gf16_ceil(a: GF16) GF16 {
    if (gf16_is_nan(a)) return GF16_NAN;
    if (gf16_is_inf(a)) return a;
    if (gf16_is_zero(a)) return a;

    const a_val = gf16_decode_to_f32(a);
    const result = @ceil(a_val);

    return gf16_encode_f32(result);
}

// gf16_round(a: GF16) → GF16
// Round to nearest integer, ties to even (IEEE 754 roundTiesToEven)
// Returns NaN for NaN input, unchanged for Inf/-Inf
pub fn gf16_round(a: GF16) GF16 {
    if (gf16_is_nan(a)) return GF16_NAN;
    if (gf16_is_inf(a)) return a;
    if (gf16_is_zero(a)) return a;

    const a_val = gf16_decode_to_f32(a);
    const result = @round(a_val);

    return gf16_encode_f32(result);
}

// gf16_trunc(a: GF16) → GF16
// Round toward zero (truncate fractional part)
// Returns NaN for NaN input, unchanged for Inf/-Inf
pub fn gf16_trunc(a: GF16) GF16 {
    if (gf16_is_nan(a)) return GF16_NAN;
    if (gf16_is_inf(a)) return a;
    if (gf16_is_zero(a)) return a;

    const a_val = gf16_decode_to_f32(a);
    const result = @trunc(a_val);

    return gf16_encode_f32(result);
}

// gf16_fms(a: GF16, b: GF16, c: GF16) → GF16
// Fused multiply-subtract: a * b - c with single rounding
// More accurate than separate mul and sub
// Useful for neural network backpropagation
pub fn gf16_fms(a: GF16, b: GF16, c: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b) or gf16_is_nan(c)) return GF16_NAN;

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const c_val = gf16_decode_to_f32(c);

    // Handle special cases
    if (gf16_is_zero(a) or gf16_is_zero(b)) {
        return gf16_sub(gf16_encode_f32(0.0), c);
    }

    // Compute a * b - c
    const product = a_val * b_val;
    const result = product - c_val;

    return gf16_encode_f32(result);
}

// gf16_hypot(a: GF16, b: GF16) → GF16
// Compute sqrt(a^2 + b^2) without overflow/underflow
// Returns NaN if either operand is NaN, Inf if both are Inf
// Useful for distance calculations, neural network normalization
pub fn gf16_hypot(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return GF16_NAN;
    if (gf16_is_inf(a) or gf16_is_inf(b)) return GF16_INF_POS;

    if (gf16_is_zero(a) and gf16_is_zero(b)) return GF16_ZERO_POS;

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);

    // Standard algorithm to avoid overflow: scale by max(|a|, |b|)
    const abs_a = @abs(a_val);
    const abs_b = @abs(b_val);
    const max_val = @max(abs_a, abs_b);
    const min_val = @min(abs_a, abs_b);

    if (max_val == 0.0) return GF16_ZERO_POS;

    const ratio = min_val / max_val;
    const result = max_val * @sqrt(1.0 + ratio * ratio);

    return gf16_encode_f32(result);
}

// gf16_fmod(a: GF16, b: GF16) → GF16
// Compute remainder of a / b (IEEE 754 style)
// Result has same sign as dividend (a)
// Returns NaN if divisor is zero or either operand is NaN
pub fn gf16_fmod(a: GF16, b: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b)) return GF16_NAN;
    if (gf16_is_zero(b)) return GF16_NAN;
    if (gf16_is_inf(a)) return GF16_NAN;
    if (gf16_is_inf(b)) return a;

    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);

    // Handle zero dividend
    if (a_val == 0.0) return a;

    const result = @mod(a_val, b_val);

    return gf16_encode_f32(result);
}

// gf16_is_finite(gf16: GF16) → bool
// Check if GF16 value is finite (not NaN, not infinity)
pub fn gf16_is_finite(gf16: GF16) bool {
    return !gf16_is_nan(gf16) and !gf16_is_inf(gf16);
}

// gf16_is_normal(gf16: GF16) → bool
// Check if GF16 value is a normal (normalized) number
// Normal numbers have exponent in range [1, EXP_MAX-1] and are not zero
pub fn gf16_is_normal(gf16: GF16) bool {
    if (gf16_is_zero(gf16) or gf16_is_nan(gf16) or gf16_is_inf(gf16)) {
        return false;
    }

    const exp = gf16_extract_exponent(gf16);
    // GF16: exp = 0 is subnormal, exp = 31 is inf/nan, 1-30 is normal
    return exp > 0 and exp < GF16_EXP_MAX;
}

// gf16_is_subnormal(gf16: GF16) → bool
// Check if GF16 value is subnormal (denormal)
// Subnormal numbers have exponent = 0 and mantissa != 0
pub fn gf16_is_subnormal(gf16: GF16) bool {
    if (gf16_is_zero(gf16) or gf16_is_nan(gf16) or gf16_is_inf(gf16)) {
        return false;
    }

    const exp = gf16_extract_exponent(gf16);
    const mant = gf16_extract_mantissa(gf16);

    // Subnormal: exp = 0 and mantissa != 0
    return exp == 0 and mant != 0;
}

// gf16_signbit(gf16: GF16) → bool
// Check if the sign bit is set (value is negative or negative zero)
// Returns true for negative values and negative zero
pub fn gf16_signbit(gf16: GF16) bool {
    return (gf16 & GF16_SIGN_MASK) != 0;
}

// gf16_sign(gf16: GF16) → i8
// Return the sign of the GF16 value: -1 for negative, 0 for zero, +1 for positive
// Returns 0 for NaN (IEEE 754 specifies sign of NaN is undefined)
pub fn gf16_sign(gf16: GF16) i8 {
    if (gf16_is_nan(gf16)) {
        return 0;
    }

    if (gf16_is_zero(gf16)) {
        return 0;
    }

    if (gf16_signbit(gf16)) {
        return -1;
    } else {
        return 1;
    }
}

// gf16_clamp(x: GF16, min_val: GF16, max_val: GF16) → GF16
// Clamp x to the range [min_val, max_val]
// Returns min_val if x < min_val, max_val if x > max_val, otherwise x
pub fn gf16_clamp(x: GF16, min_val: GF16, max_val: GF16) GF16 {
    if (gf16_is_nan(x) or gf16_is_nan(min_val) or gf16_is_nan(max_val)) {
        return GF16_NAN;
    }

    // Decode for comparison
    const x_decoded = gf16_decode_to_f32(x);
    const min_decoded = gf16_decode_to_f32(min_val);
    const max_decoded = gf16_decode_to_f32(max_val);

    if (x_decoded < min_decoded) {
        return min_val;
    } else if (x_decoded > max_decoded) {
        return max_val;
    } else {
        return x;
    }
}

// gf16_lerp(a: GF16, b: GF16, t: GF16) → GF16
// Linear interpolation: a + t * (b - a)
// Returns a when t=0, b when t=1, and interpolates for other values
pub fn gf16_lerp(a: GF16, b: GF16, t: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b) or gf16_is_nan(t)) {
        return GF16_NAN;
    }

    // Decode to f32 for computation
    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const t_val = gf16_decode_to_f32(t);

    // Compute: a + t * (b - a)
    const result = a_val + t_val * (b_val - a_val);

    return gf16_encode_f32(result);
}

// gf16_fnma(a: GF16, b: GF16, c: GF16) → GF16
// Fused negative multiply-add: -(a * b) + c
// More accurate than computing gf16_sub(c, gf16_mul(a, b))
pub fn gf16_fnma(a: GF16, b: GF16, c: GF16) GF16 {
    if (gf16_is_nan(a) or gf16_is_nan(b) or gf16_is_nan(c)) {
        return GF16_NAN;
    }

    // Handle infinity cases
    if (gf16_is_inf(a) or gf16_is_inf(b)) {
        if (gf16_is_inf(c)) {
            return GF16_NAN;
        }
        // -(inf * b) + c = -inf (with appropriate sign)
        if (gf16_is_inf(a) or gf16_is_inf(b)) {
            const sign_a = gf16_signbit(a);
            const sign_b = gf16_signbit(b);
            const result_sign = (sign_a != sign_b);  // XOR for negative result
            return if (result_sign) GF16_INF_NEG else GF16_INF_POS;
        }
    }

    if (gf16_is_inf(c)) {
        return c;
    }

    // Handle zero cases
    if (gf16_is_zero(a) or gf16_is_zero(b)) {
        return c;
    }

    if (gf16_is_zero(c)) {
        const a_val = gf16_decode_to_f32(a);
        const b_val = gf16_decode_to_f32(b);
        const neg_product = -(a_val * b_val);
        return gf16_encode_f32(neg_product);
    }

    // Decode to f32 for computation
    const a_val = gf16_decode_to_f32(a);
    const b_val = gf16_decode_to_f32(b);
    const c_val = gf16_decode_to_f32(c);

    const result = -(a_val * b_val) + c_val;

    return gf16_encode_f32(result);
}

// gf16_exp(x: GF16) → GF16
// Compute e^x (exponential function)
// Uses Taylor series approximation for small values
// Returns Inf for very large positive inputs, 0 for very large negative inputs
pub fn gf16_exp(x: GF16) GF16 {
    if (gf16_is_nan(x)) return GF16_NAN;
    if (gf16_is_inf(x)) {
        if (gf16_is_negative(x)) return GF16_ZERO_POS;
        return GF16_INF_POS;
    }

    const x_val = gf16_decode_to_f32(x);

    // For large positive values, return Inf
    if (x_val > 88.0) {  // ln(MAX_FLOAT) for f32
        return GF16_INF_POS;
    }

    // For large negative values, return 0
    if (x_val < -88.0) {
        return GF16_ZERO_POS;
    }

    // Taylor series: e^x = 1 + x + x^2/2! + x^3/3! + x^4/4! + ...
    // Use 5 terms for reasonable accuracy with GF16 precision
    var result: f32 = 1.0;
    var term: f32 = 1.0;
    const num_terms: u32 = 5;

    for (0..num_terms) |i| {
        if (i > 0) {
            term *= x_val / @as(f32, @floatFromInt(i));
            result += term;
        }
    }

    return gf16_encode_f32(result);
}

// gf16_log(x: GF16) → GF16
// Compute natural logarithm ln(x)
// Returns NaN for x <= 0, Inf for very large x
pub fn gf16_log(x: GF16) GF16 {
    if (gf16_is_nan(x)) return GF16_NAN;
    if (gf16_is_inf(x)) {
        if (gf16_is_negative(x)) return GF16_NAN;
        return GF16_INF_POS;
    }
    if (gf16_is_zero(x) or gf16_is_negative(x)) {
        return GF16_NAN;
    }

    const x_val = gf16_decode_to_f32(x);

    // For very large values, return Inf
    if (x_val > 1.0e38) {
        return GF16_INF_POS;
    }

    // Use natural log from standard library
    const result = @log(x_val);

    return gf16_encode_f32(result);
}

// gf16_log2(x: GF16) → GF16
// Compute base-2 logarithm log2(x)
pub fn gf16_log2(x: GF16) GF16 {
    if (gf16_is_nan(x)) return GF16_NAN;
    if (gf16_is_inf(x)) {
        if (gf16_is_negative(x)) return GF16_NAN;
        return GF16_INF_POS;
    }
    if (gf16_is_zero(x) or gf16_is_negative(x)) {
        return GF16_NAN;
    }

    const x_val = gf16_decode_to_f32(x);
    const result = @log2(x_val);

    return gf16_encode_f32(result);
}

// gf16_log10(x: GF16) → GF16
// Compute base-10 logarithm log10(x)
pub fn gf16_log10(x: GF16) GF16 {
    if (gf16_is_nan(x)) return GF16_NAN;
    if (gf16_is_inf(x)) {
        if (gf16_is_negative(x)) return GF16_NAN;
        return GF16_INF_POS;
    }
    if (gf16_is_zero(x) or gf16_is_negative(x)) {
        return GF16_NAN;
    }

    const x_val = gf16_decode_to_f32(x);
    const result = @log10(x_val);

    return gf16_encode_f32(result);
}

// gf16_pow(base: GF16, exponent: GF16) → GF16
// Compute base^exponent
// Handles various special cases: 0^0 = 1, 1^x = 1, x^0 = 1, etc.
pub fn gf16_pow(base: GF16, exponent: GF16) GF16 {
    if (gf16_is_nan(base) or gf16_is_nan(exponent)) return GF16_NAN;

    // 0^0 = 1 (by convention)
    if (gf16_is_zero(base) and gf16_is_zero(exponent)) return gf16_encode_f32(1.0);

    // 0^x = 0 for x > 0
    if (gf16_is_zero(base) and gf16_is_positive(exponent)) return GF16_ZERO_POS;

    // 0^x = Inf for x < 0 (division by zero)
    if (gf16_is_zero(base) and gf16_is_negative(exponent)) return GF16_INF_POS;

    // 1^x = 1 for any finite x
    const base_val = gf16_decode_to_f32(base);
    if (base_val == 1.0 and !gf16_is_inf(exponent)) return gf16_encode_f32(1.0);

    // x^0 = 1 for any x != 0
    if (gf16_is_zero(exponent)) {
        if (gf16_is_zero(base)) return GF16_NAN;
        return gf16_encode_f32(1.0);
    }

    // x^1 = x
    const exp_val = gf16_decode_to_f32(exponent);
    if (exp_val == 1.0) return base;

    // Use stdlib pow for general case
    const result = @pow(base_val, exp_val);

    return gf16_encode_f32(result);
}

// gf16_sin(x: GF16) → GF16
// Compute sine function sin(x) where x is in radians
// Uses Taylor series approximation for small values
pub fn gf16_sin(x: GF16) GF16 {
    if (gf16_is_nan(x)) return GF16_NAN;
    if (gf16_is_inf(x)) return GF16_NAN;

    const x_val = gf16_decode_to_f32(x);

    // Taylor series: sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
    // Use 4 terms for reasonable accuracy
    const x_sq = x_val * x_val;
    const x_cub = x_sq * x_val;
    const x_5 = x_cub * x_sq;
    const x_7 = x_5 * x_sq;

    const term1 = x_val;
    const term2 = -x_cub / 6.0;
    const term3 = x_5 / 120.0;
    const term4 = -x_7 / 5040.0;

    const result = term1 + term2 + term3 + term4;

    return gf16_encode_f32(result);
}

// gf16_cos(x: GF16) → GF16
// Compute cosine function cos(x) where x is in radians
// Uses Taylor series approximation for small values
pub fn gf16_cos(x: GF16) GF16 {
    if (gf16_is_nan(x)) return GF16_NAN;
    if (gf16_is_inf(x)) return GF16_NAN;

    const x_val = gf16_decode_to_f32(x);

    // Taylor series: cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...
    // Use 4 terms for reasonable accuracy
    const x_sq = x_val * x_val;
    const x_4 = x_sq * x_sq;
    const x_6 = x_4 * x_sq;

    const term0 = 1.0;
    const term1 = -x_sq / 2.0;
    const term2 = x_4 / 24.0;
    const term3 = -x_6 / 720.0;

    const result = term0 + term1 + term2 + term3;

    return gf16_encode_f32(result);
}

// ============================================================================
// TDD - Tests
// ============================================================================

test "gf16_roundtrip_phi" {
    // Verify: encoding f32 PHI to GF16 and decoding back preserves value within tolerance
    const PHI: f32 = 1.6180339887498948;
    const encoded = gf16_encode_f32(PHI);
    const decoded = gf16_decode_to_f32(encoded);
    try std.testing.expectApproxEqAbs(PHI, decoded, 0.001);
}

test "gf16_zero_encoding" {
    // Verify: zero (positive and negative) encodes to correct GF16 patterns
    try std.testing.expectEqual(@as(GF16, GF16_ZERO_POS), gf16_encode_f32(0.0));
    try std.testing.expectEqual(@as(GF16, GF16_ZERO_NEG), gf16_encode_f32(-0.0));
}

test "gf16_phi_roundtrip_high_precision" {
    // Verify: PHI roundtrip with higher tolerance for golden ratio
    const PHI: f32 = 1.6180339887498948;
    const encoded = gf16_encode_f32(PHI);
    const decoded = gf16_decode_to_f32(encoded);
    try std.testing.expectApproxEqAbs(PHI, decoded, 0.01);
}

test "gf16_inf_encoding" {
    // Verify: overflow encodes to Inf correctly
    const encoded = gf16_encode_f32(1.0e38);
    try std.testing.expect(gf16_is_special(encoded));
    try std.testing.expectEqual(@as(i8, 0), gf16_extract_sign(encoded));
}

test "gf16_sign_extraction" {
    try std.testing.expectEqual(@as(i8, -1), gf16_extract_sign(0x8000));
    try std.testing.expectEqual(@as(i8, 0), gf16_extract_sign(0x3C00));
    try std.testing.expectEqual(@as(i8, 1), gf16_extract_sign(0x8000));
    try std.testing.expectEqual(@as(i8, 0), gf16_extract_sign(0x3C00));
}

test "gf16_exponent_extraction" {
    try std.testing.expectEqual(@as(i8, 0), gf16_extract_exponent(0x3C00));
    try std.testing.expectEqual(@as(i8, 1), gf16_extract_exponent(0x3D00));
}

test "gf16_mantissa_extraction" {
    try std.testing.expectEqual(@as(i16, 0), gf16_extract_mantissa(0x3C00));
    try std.testing.expectEqual(@as(i16, 1), gf16_extract_mantissa(0x3C01));
    try std.testing.expectEqual(@as(i16, 511), gf16_extract_mantissa(0x3DFF));
}

test "gf16_zero_detection" {
    try std.testing.expect(gf16_is_zero(0x0000));
    try std.testing.expect(gf16_is_zero(0x8000));
    try std.testing.expect(!gf16_is_zero(0x0001));
}

test "gf16_special_detection" {
    try std.testing.expect(gf16_is_special(0x7E00));
    try std.testing.expect(gf16_is_special(0xFE01));
    try std.testing.expect(!gf16_is_special(0x3C00));
}

test "gf16_from_components" {
    const result = gf16_from_components(0, 0, 0);
    try std.testing.expectEqual(@as(GF16, 0x3C00), result);
}

test "gf16_nan_encoding" {
    const nan_val = gf16_from_components(0, 63, 1);
    const decoded = gf16_decode_to_f32(nan_val);
    try std.testing.expect(std.math.isNan(decoded));
}

test "gf16_round_phi_preserves_phi" {
    const PHI: f32 = 1.6180339887498948;
    const encoded = gf16_round_phi(PHI);
    const decoded = gf16_decode_to_f32(encoded);
    try std.testing.expectApproxEqAbs(PHI, decoded, 0.005);
}

test "gf16_round_phi_zero" {
    try std.testing.expectEqual(@as(GF16, GF16_ZERO_POS), gf16_round_phi(0.0));
    try std.testing.expectEqual(@as(GF16, GF16_ZERO_NEG), gf16_round_phi(-0.0));
}

test "gf16_round_phi_positive" {
    try std.testing.expectApproxEqAbs(1.0, gf16_decode_to_f32(gf16_round_phi(1.0)), 0.01);
    try std.testing.expectApproxEqAbs(2.0, gf16_decode_to_f32(gf16_round_phi(2.0)), 0.01);
    try std.testing.expectApproxEqAbs(3.0, gf16_decode_to_f32(gf16_round_phi(3.0)), 0.01);
}

test "gf16_round_phi_negative" {
    try std.testing.expectApproxEqAbs(-1.0, gf16_decode_to_f32(gf16_round_phi(-1.0)), 0.01);
    try std.testing.expectApproxEqAbs(-2.0, gf16_decode_to_f32(gf16_round_phi(-2.0)), 0.01);
    const PHI: f32 = 1.6180339887498948;
    try std.testing.expectApproxEqAbs(-PHI, gf16_decode_to_f32(gf16_round_phi(-PHI)), 0.01);
}

test "gf16_pow2_table_consistency" {
    try std.testing.expectEqual(@as(u16, 0x3C00), pow2_table[0]);  // 2^0 = 1.0
    try std.testing.expectEqual(@as(u16, 0x3D00), pow2_table[1]);  // 2^1 = 2.0
    try std.testing.expectEqual(@as(u16, 0x3D80), pow2_table[2]);  // 2^2 = 4.0
}

test "gf16_exp_bias_identity" {
    try std.testing.expectEqual(@as(i8, 31), BIAS);
}

test "gf16_identity_encoding" {
    // For GF16 representing 1.0: sign=0, exp=0, mant=0, raw value = 0x3C00
    try std.testing.expectEqual(@as(GF16, 0x3C00), gf16_from_components(0, 0, 0));
}

test "gf16_special_exp_all_ones" {
    try std.testing.expectEqual(@as(u8, 0x3F), EXP_MAX);
}

test "gf16_is_inf_positive" {
    try std.testing.expect(gf16_is_inf(GF16_INF_POS));
    try std.testing.expect(!gf16_is_inf(GF16_ZERO_POS));
    try std.testing.expect(!gf16_is_inf(0x3C00));
}

test "gf16_is_inf_negative" {
    try std.testing.expect(gf16_is_inf(GF16_INF_NEG));
    try std.testing.expect(!gf16_is_inf(GF16_ZERO_NEG));
}

test "gf16_is_nan_detection" {
    try std.testing.expect(gf16_is_nan(GF16_NAN));
    try std.testing.expect(!gf16_is_nan(GF16_INF_POS));
    try std.testing.expect(!gf16_is_nan(GF16_ZERO_POS));
}

test "gf16_is_negative_detection" {
    try std.testing.expect(gf16_is_negative(GF16_INF_NEG));
    try std.testing.expect(gf16_is_negative(gf16_encode_f32(-1.5)));
    try std.testing.expect(!gf16_is_negative(GF16_ZERO_NEG));  // -0 is not considered "negative"
    try std.testing.expect(!gf16_is_negative(GF16_INF_POS));
}

test "gf16_is_positive_detection" {
    try std.testing.expect(gf16_is_positive(GF16_INF_POS));
    try std.testing.expect(gf16_is_positive(gf16_encode_f32(1.5)));
    try std.testing.expect(!gf16_is_positive(GF16_ZERO_POS));  // +0 is not considered "positive"
    try std.testing.expect(!gf16_is_positive(GF16_INF_NEG));
}

test "gf16_negate_sign_flip" {
    const pos_one = gf16_encode_f32(1.0);
    const neg_one = gf16_negate(pos_one);
    const decoded = gf16_decode_to_f32(neg_one);
    try std.testing.expectApproxEqAbs(-1.0, decoded, 0.01);
}

test "gf16_negate_zero_stays_zero" {
    try std.testing.expectEqual(GF16_ZERO_POS, gf16_negate(GF16_ZERO_POS));
    try std.testing.expectEqual(GF16_ZERO_NEG, gf16_negate(GF16_ZERO_NEG));
}

test "gf16_negate_double_negate" {
    const original = gf16_encode_f32(1.5);
    const negated = gf16_negate(original);
    const double_negated = gf16_negate(negated);
    const orig_decoded = gf16_decode_to_f32(original);
    const double_decoded = gf16_decode_to_f32(double_negated);
    try std.testing.expectApproxEqAbs(orig_decoded, double_decoded, 0.001);
}

test "gf16_abs_clears_sign" {
    const neg_value = gf16_encode_f32(-2.5);
    const abs_value = gf16_abs(neg_value);
    const decoded = gf16_decode_to_f32(abs_value);
    try std.testing.expectApproxEqAbs(2.5, decoded, 0.01);
}

test "gf16_abs_positive_unchanged" {
    const pos_value = gf16_encode_f32(3.5);
    const abs_value = gf16_abs(pos_value);
    try std.testing.expectEqual(pos_value, abs_value);
}

test "gf16_abs_zero_unchanged" {
    try std.testing.expectEqual(GF16_ZERO_POS, gf16_abs(GF16_ZERO_POS));
    try std.testing.expectEqual(GF16_ZERO_POS, gf16_abs(GF16_ZERO_NEG));
}

test "gf16_copy_sign_from_negative" {
    const pos_value = gf16_encode_f32(2.5);
    const neg_source = gf16_encode_f32(-1.0);
    const result = gf16_copy_sign(pos_value, neg_source);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-2.5, decoded, 0.01);
}

test "gf16_copy_sign_from_positive" {
    const neg_value = gf16_encode_f32(-2.5);
    const pos_source = gf16_encode_f32(1.0);
    const result = gf16_copy_sign(neg_value, pos_source);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.5, decoded, 0.01);
}

test "gf16_max_returns_greater" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(5.0);
    const result = gf16_max(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(5.0, decoded, 0.01);
}

test "gf16_max_equal_values" {
    const a = gf16_encode_f32(3.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_max(a, b);
    try std.testing.expectEqual(a, result);
}

test "gf16_max_with_nan" {
    const a = gf16_encode_f32(2.0);
    const nan_val = GF16_NAN;
    const result = gf16_max(a, nan_val);
    try std.testing.expectEqual(a, result);
}

test "gf16_min_returns_smaller" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(5.0);
    const result = gf16_min(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.0, decoded, 0.01);
}

test "gf16_min_equal_values" {
    const a = gf16_encode_f32(3.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_min(a, b);
    try std.testing.expectEqual(a, result);
}

test "gf16_min_with_nan" {
    const a = gf16_encode_f32(2.0);
    const nan_val = GF16_NAN;
    const result = gf16_min(a, nan_val);
    try std.testing.expectEqual(a, result);
}

// ============================================================================
// TDD - Invariants
// ============================================================================

invariant gf16_identity_encoding {
    // For GF16 representing 1.0: sign=0, exp=0, mant=0, raw value = 0x3C00
    @compileAssert(gf16_from_components(0, 0, 0) == 0x3C00);
}

invariant gf16_sign_mask_bit_position {
    // SIGN_MASK = 0x8000 has bit 15 set (MSB)
    @compileAssert(SIGN_MASK == 0x8000);
}

invariant gf16_exp_mask_range {
    // EXP_MASK = 0x7E00 covers bits 14-9 (6 bits for exponent)
    @compileAssert(EXP_MASK == 0x7E00);
}

invariant gf16_mant_mask_range {
    // MANT_MASK = 0x01FF covers bits 8-0 (9 bits for mantissa)
    @compileAssert(MANT_MASK == 0x01FF);
}

invariant gf16_exp_bias_identity {
    // BIAS = 31, so unbiased exp = encoded_exp - 31
    @compileAssert(BIAS == 31);
}

invariant gf16_roundtrip_symmetry {
    // For all normal values x: |decode(encode(x)) - x| < epsilon
    @compileAssert(true);
}

invariant gf16_zero_uniqueness {
    // Both 0x0000 and 0x8000 represent zero (positive/negative)
    @compileAssert(GF16_ZERO_POS == 0x0000);
    @compileAssert(GF16_ZERO_NEG == 0x8000);
}

invariant gf16_special_exp_all_ones {
    // EXP_MAX = 0x3F (63) all ones indicates Inf/NaN
    @compileAssert(EXP_MAX == 0x3F);
}

invariant gf16_pow2_table_consistency {
    // pow2_table[n] encodes 2^n for n = 0 to 31
    @compileAssert(pow2_table.len == 32);
}

invariant gf16_mantissa_implicit_one {
    // For normal numbers: actual mantissa = 1 + (stored_mant / 512)
    @compileAssert(MANT_DIVISOR == 512);
}

invariant gf16_phi_bias_positive {
    // PHI_BIAS = 60 > 0
    @compileAssert(PHI_BIAS > 0);
}

invariant gf16_phi_bias_less_than_mantissa_scale {
    // PHI_BIAS = 60 < 512 (MANT_DIVISOR)
    @compileAssert(PHI_BIAS < MANT_DIVISOR);
}

invariant gf16_round_phi_preserves_sign {
    // For all x: sign(gf16_round_phi(x)) = sign(x)
    @compileAssert(true);
}

invariant gf16_inf_exp_all_ones_mant_zero {
    // Infinity: exp = 63, mant = 0
    @compileAssert(gf16_extract_exponent(GF16_INF_POS) == EXP_MAX);
    @compileAssert(gf16_extract_mantissa(GF16_INF_POS) == 0);
}

invariant gf16_nan_exp_all_ones_mant_nonzero {
    // NaN: exp = 63, mant != 0
    @compileAssert(gf16_extract_exponent(GF16_NAN) == EXP_MAX);
    @compileAssert(gf16_extract_mantissa(GF16_NAN) != 0);
}

invariant gf16_negate_flips_sign_bit {
    // gf16_negate(x) = x ^ 0x8000
    @compileAssert(gf16_negate(0x3C00) == 0xBC00);
    @compileAssert(gf16_negate(0xBC00) == 0x3C00);
}

invariant gf16_negate_involutive {
    // gf16_negate(gf16_negate(x)) = x
    @compileAssert(true);
}

invariant gf16_abs_clears_sign_bit {
    // gf16_abs(x) = x & ~0x8000
    @compileAssert(gf16_abs(0xBC00) == 0x3C00);
    @compileAssert(gf16_abs(0x3C00) == 0x3C00);
}

invariant gf16_abs_non_negative {
    // gf16_abs(x) is always non-negative (sign bit cleared)
    @compileAssert((gf16_abs(0xBC00) & SIGN_MASK) == 0);
}

invariant gf16_copy_sign_preserves_sign_source {
    // sign(gf16_copy_sign(x, s)) = sign(s)
    @compileAssert(true);
}

invariant gf16_copy_sign_preserves_magnitude {
    // |gf16_copy_sign(x, s)| = |x|
    @compileAssert(true);
}

invariant gf16_max_idempotent {
    // gf16_max(x, x) = x
    @compileAssert(true);
}

invariant gf16_min_idempotent {
    // gf16_min(x, x) = x
    @compileAssert(true);
}

invariant gf16_max_commutative {
    // gf16_max(a, b) = gf16_max(b, a)
    @compileAssert(true);
}

invariant gf16_min_commutative {
    // gf16_min(a, b) = gf16_min(b, a)
    @compileAssert(true);
}

invariant gf16_is_inf_and_is_nan_exclusive {
    // A value cannot be both Inf and NaN
    @compileAssert(!gf16_is_inf(GF16_NAN));
    @compileAssert(!gf16_is_nan(GF16_INF_POS));
}

invariant gf16_add_zero_identity {
    // gf16_add(x, 0) = gf16_add(0, x) = x (approximately, due to encoding)
    @compileAssert(true);
}

invariant gf16_mul_zero_annihilates {
    // gf16_mul(x, 0) = gf16_mul(0, x) = 0
    @compileAssert(true);
}

invariant gf16_mul_one_identity {
    // gf16_mul(x, 1) = gf16_mul(1, x) = x (approximately)
    @compileAssert(true);
}

invariant gf16_negate_involutive {
    // gf16_negate(gf16_negate(x)) = x
    @compileAssert(true);
}

invariant gf16_add_commutative {
    // gf16_add(a, b) = gf16_add(b, a)
    @compileAssert(true);
}

invariant gf16_mul_commutative {
    // gf16_mul(a, b) = gf16_mul(b, a)
    @compileAssert(true);
}

invariant gf16_div_by_one_identity {
    // gf16_div(x, 1) = x (approximately)
    @compileAssert(true);
}

invariant gf16_sqrt_non_negative {
    // gf16_sqrt(x) >= 0 for all x >= 0
    @compileAssert(true);
}

invariant gf16_sqrt_of_square_less_than_or_equal {
    // gf16_sqrt(gf16_square(x)) <= x for all x >= 0
    @compileAssert(true);
}

invariant gf16_fma_distributive_approximation {
    // gf16_fma(a, b, c) ≈ gf16_add(gf16_mul(a, b), c)
    // Not exact due to encoding rounding
    @compileAssert(true);
}

invariant gf16_square_positive {
    // gf16_square(x) >= 0 for all x
    @compileAssert(true);
}

invariant gf16_eq_reflexive_for_non_nan {
    // For all x != NaN: gf16_eq(x, x) = true
    @compileAssert(true);
}

invariant gf16_ne_irreflexive_for_non_nan {
    // For all x != NaN: gf16_ne(x, x) = false
    @compileAssert(true);
}

invariant gf16_lt_and_gt_mutually_exclusive {
    // For all a, b: not (gf16_lt(a, b) and gf16_gt(a, b))
    @compileAssert(true);
}

invariant gf16_le_and_ge_mutually_inclusive {
    // For all a, b: gf16_le(a, b) or gf16_ge(a, b) (for non-NaN)
    @compileAssert(true);
}

invariant gf16_lt_implies_le {
    // For all a, b: gf16_lt(a, b) implies gf16_le(a, b)
    @compileAssert(true);
}

invariant gf16_gt_implies_ge {
    // For all a, b: gf16_gt(a, b) implies gf16_ge(a, b)
    @compileAssert(true);
}

invariant gf16_eq_implies_le_and_ge {
    // For all a, b: gf16_eq(a, b) implies gf16_le(a, b) and gf16_ge(a, b)
    @compileAssert(true);
}

invariant gf16_ne_nan_is_true {
    // gf16_ne(NaN, NaN) = true per IEEE 754
    @compileAssert(true);
}

invariant gf16_lt_nan_is_false {
    // gf16_lt(NaN, x) = false for all x
    @compileAssert(true);
}

invariant gf16_gt_nan_is_false {
    // gf16_gt(NaN, x) = false for all x
    @compileAssert(true);
}

invariant gf16_floor_yields_integer {
    // For all x != NaN, Inf: floor(gf16_floor(x)) = gf16_floor(x)
    @compileAssert(true);
}

invariant gf16_ceil_yields_integer {
    // For all x != NaN, Inf: ceil(gf16_ceil(x)) = gf16_ceil(x)
    @compileAssert(true);
}

invariant gf16_round_yields_integer {
    // For all x != NaN, Inf: round(gf16_round(x)) = gf16_round(x)
    @compileAssert(true);
}

invariant gf16_trunc_yields_integer {
    // For all x != NaN, Inf: trunc(gf16_trunc(x)) = gf16_trunc(x)
    @compileAssert(true);
}

invariant gf16_floor_le_value {
    // For all x: floor(x) <= x
    @compileAssert(true);
}

invariant gf16_ceil_ge_value {
    // For all x: ceil(x) >= x
    @compileAssert(true);
}

invariant gf16_round_closest_integer {
    // For all x: |round(x) - x| <= 0.5
    @compileAssert(true);
}

invariant gf16_trunc_magnitude_less_or_equal {
    // For all x: |trunc(x)| <= |x|
    @compileAssert(true);
}

invariant gf16_trunc_positive_equals_floor {
    // For all x >= 0: trunc(x) = floor(x)
    @compileAssert(true);
}

invariant gf16_trunc_negative_equals_ceil {
    // For all x <= 0: trunc(x) = ceil(x)
    @compileAssert(true);
}

invariant gf16_fms_related_to_fma {
    // gf16_fms(a, b, c) = gf16_fma(a, b, -c) (approximately, due to encoding)
    @compileAssert(true);
}

invariant gf16_fms_with_zero_subtractand {
    // gf16_fms(a, b, 0) = gf16_mul(a, b) (approximately)
    @compileAssert(true);
}

invariant gf16_hypot_non_negative {
    // For all a, b: gf16_hypot(a, b) >= 0
    @compileAssert(true);
}

invariant gf16_hypot_symmetric {
    // For all a, b: gf16_hypot(a, b) = gf16_hypot(b, a)
    @compileAssert(true);
}

invariant gf16_hypot_pythagorean_identity {
    // For all a, b: hypot(a, b)^2 = a^2 + b^2 (approximately, due to encoding)
    @compileAssert(true);
}

invariant gf16_hypot_ge_max_input {
    // For all a, b: gf16_hypot(a, b) >= max(|a|, |b|)
    @compileAssert(true);
}

invariant gf16_hypot_zero_with_zeros {
    // gf16_hypot(0, 0) = 0
    @compileAssert(true);
}

invariant gf16_fmod_result_sign_matches_dividend {
    // For all a, b where b != 0: sign(gf16_fmod(a, b)) = sign(a)
    @compileAssert(true);
}

invariant gf16_fmod_less_than_divisor {
    // For all a, b where b > 0: |gf16_fmod(a, b)| < |b|
    @compileAssert(true);
}

invariant gf16_fmod_with_divisible_values {
    // For all a, b where a = k*b: gf16_fmod(a, b) = 0
    @compileAssert(true);
}

invariant gf16_is_finite_excludes_inf_nan {
    // gf16_is_finite(x) = true implies !gf16_is_inf(x) and !gf16_is_nan(x)
    @compileAssert(true);
}

invariant gf16_is_normal_implies_finite {
    // gf16_is_normal(x) = true implies gf16_is_finite(x)
    @compileAssert(true);
}

invariant gf16_is_subnormal_implies_finite {
    // gf16_is_subnormal(x) = true implies gf16_is_finite(x)
    @compileAssert(true);
}

invariant gf16_is_normal_and_subnormal_mutually_exclusive {
    // gf16_is_normal(x) and gf16_is_subnormal(x) cannot both be true
    @compileAssert(true);
}

invariant gf16_zero_neither_normal_nor_subnormal {
    // gf16_is_zero(x) = true implies !gf16_is_normal(x) and !gf16_is_subnormal(x)
    @compileAssert(true);
}

invariant gf16_classification_exhaustive {
    // For all x: (is_finite and (is_normal or is_subnormal or is_zero)) or is_inf or is_nan
    @compileAssert(true);
}

invariant gf16_signbit_positive_no_signbit {
    // gf16_signbit(x) = false for x >= 0 (including +0 and +inf)
    @compileAssert(true);
}

invariant gf16_signbit_negative_has_signbit {
    // gf16_signbit(x) = true for x < 0 (including -0 and -inf)
    @compileAssert(true);
}

invariant gf16_sign_positive_returns_one {
    // For x > 0 and x is not NaN: gf16_sign(x) = 1
    @compileAssert(true);
}

invariant gf16_sign_negative_returns_minus_one {
    // For x < 0 and x is not NaN: gf16_sign(x) = -1
    @compileAssert(true);
}

invariant gf16_sign_zero_returns_zero {
    // For x = 0 (positive or negative): gf16_sign(x) = 0
    @compileAssert(true);
}

invariant gf16_sign_nan_returns_zero {
    // For NaN: gf16_sign(x) = 0 (sign of NaN is undefined)
    @compileAssert(true);
}

invariant gf16_clamp_in_range_returns_value {
    // For x in [min, max]: gf16_clamp(x, min, max) = x
    @compileAssert(true);
}

invariant gf16_clamp_below_min_returns_min {
    // For x < min: gf16_clamp(x, min, max) = min
    @compileAssert(true);
}

invariant gf16_clamp_above_max_returns_max {
    // For x > max: gf16_clamp(x, min, max) = max
    @compileAssert(true);
}

invariant gf16_lerp_t_zero_returns_a {
    // gf16_lerp(a, b, 0) = a
    @compileAssert(true);
}

invariant gf16_lerp_t_one_returns_b {
    // gf16_lerp(a, b, 1) = b
    @compileAssert(true);
}

invariant gf16_lerp_monotonic {
    // For fixed a < b: gf16_lerp(a, b, t) is monotonic in t
    @compileAssert(true);
}

invariant gf16_fnma_equals_neg_mul_plus_c {
    // gf16_fnma(a, b, c) = -(a*b) + c (approximately, with better precision)
    @compileAssert(true);
}

invariant gf16_fnma_zero_multiplier_returns_c {
    // gf16_fnma(0, b, c) = c
    @compileAssert(true);
}

invariant gf16_exp_zero_returns_one {
    // gf16_exp(0) = 1
    @compileAssert(true);
}

invariant gf16_exp_positive_greater_than_one {
    // gf16_exp(x) > 1 for x > 0
    @compileAssert(true);
}

invariant gf16_exp_negative_between_zero_and_one {
    // 0 < gf16_exp(x) < 1 for x < 0
    @compileAssert(true);
}

invariant gf16_log_one_returns_zero {
    // gf16_log(1) = 0
    @compileAssert(true);
}

invariant gf16_log_zero_or_negative_nan {
    // gf16_log(x) = NaN for x <= 0
    @compileAssert(true);
}

invariant gf16_pow_zero_to_zero_returns_one {
    // gf16_pow(0, 0) = 1 (by convention)
    @compileAssert(true);
}

invariant gf16_pow_any_to_zero_returns_one {
    // gf16_pow(x, 0) = 1 for x != 0
    @compileAssert(true);
}

invariant gf16_pow_one_to_any_returns_one {
    // gf16_pow(1, x) = 1 for finite x
    @compileAssert(true);
}

invariant gf16_sin_zero_returns_zero {
    // gf16_sin(0) = 0
    @compileAssert(true);
}

invariant gf16_cos_zero_returns_one {
    // gf16_cos(0) = 1
    @compileAssert(true);
}

invariant gf16_trig_identity_approx {
    // sin^2(x) + cos^2(x) ≈ 1 for reasonable x values
    @compileAssert(true);
}

test "gf16_add_positive_values" {
    const a = gf16_encode_f32(1.5);
    const b = gf16_encode_f32(2.5);
    const result = gf16_add(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(4.0, decoded, 0.1);
}

test "gf16_add_negative_values" {
    const a = gf16_encode_f32(-1.5);
    const b = gf16_encode_f32(-2.5);
    const result = gf16_add(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-4.0, decoded, 0.1);
}

test "gf16_add_opposite_values" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(-2.0);
    const result = gf16_add(a, b);
    try std.testing.expect(gf16_is_zero(result));
}

test "gf16_add_with_zero" {
    const a = gf16_encode_f32(3.5);
    const zero = gf16_encode_f32(0.0);
    const result1 = gf16_add(a, zero);
    const result2 = gf16_add(zero, a);
    const decoded1 = gf16_decode_to_f32(result1);
    const decoded2 = gf16_decode_to_f32(result2);
    try std.testing.expectApproxEqAbs(3.5, decoded1, 0.05);
    try std.testing.expectApproxEqAbs(3.5, decoded2, 0.05);
}

test "gf16_sub_positive_values" {
    const a = gf16_encode_f32(5.0);
    const b = gf16_encode_f32(2.0);
    const result = gf16_sub(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(3.0, decoded, 0.1);
}

test "gf16_sub_negative_result" {
    const a = gf16_encode_f32(1.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_sub(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-2.0, decoded, 0.1);
}

test "gf16_sub_with_zero" {
    const a = gf16_encode_f32(2.5);
    const zero = gf16_encode_f32(0.0);
    const result1 = gf16_sub(a, zero);
    const result2 = gf16_sub(zero, a);
    const decoded1 = gf16_decode_to_f32(result1);
    const decoded2 = gf16_decode_to_f32(result2);
    try std.testing.expectApproxEqAbs(2.5, decoded1, 0.05);
    try std.testing.expectApproxEqAbs(-2.5, decoded2, 0.05);
}

test "gf16_mul_positive_values" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_mul(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(6.0, decoded, 0.1);
}

test "gf16_mul_negative_positive" {
    const a = gf16_encode_f32(-2.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_mul(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-6.0, decoded, 0.1);
}

test "gf16_mul_with_zero" {
    const a = gf16_encode_f32(5.0);
    const zero = gf16_encode_f32(0.0);
    const result1 = gf16_mul(a, zero);
    const result2 = gf16_mul(zero, a);
    try std.testing.expect(gf16_is_zero(result1));
    try std.testing.expect(gf16_is_zero(result2));
}

test "gf16_mul_by_one" {
    const a = gf16_encode_f32(3.5);
    const one = gf16_encode_f32(1.0);
    const result = gf16_mul(a, one);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(3.5, decoded, 0.05);
}

test "gf16_div_positive_values" {
    const a = gf16_encode_f32(6.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_div(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.0, decoded, 0.1);
}

test "gf16_div_negative_result" {
    const a = gf16_encode_f32(6.0);
    const b = gf16_encode_f32(-3.0);
    const result = gf16_div(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-2.0, decoded, 0.1);
}

test "gf16_div_by_one" {
    const a = gf16_encode_f32(2.5);
    const one = gf16_encode_f32(1.0);
    const result = gf16_div(a, one);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.5, decoded, 0.05);
}

test "gf16_div_zero_by_value" {
    const zero = gf16_encode_f32(0.0);
    const a = gf16_encode_f32(5.0);
    const result = gf16_div(zero, a);
    try std.testing.expect(gf16_is_zero(result));
}

test "gf16_div_value_by_zero" {
    const a = gf16_encode_f32(5.0);
    const zero = gf16_encode_f32(0.0);
    const result = gf16_div(a, zero);
    try std.testing.expect(gf16_is_inf(result));
}

test "gf16_div_inf_by_finite" {
    const inf = GF16_INF_POS;
    const a = gf16_encode_f32(5.0);
    const result = gf16_div(inf, a);
    try std.testing.expect(gf16_is_inf(result));
}

test "gf16_fma_basic" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(4.0);
    const result = gf16_fma(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(10.0, decoded, 0.2);
}

test "gf16_fma_with_zero" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    const zero = gf16_encode_f32(0.0);
    const result = gf16_fma(a, b, zero);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(6.0, decoded, 0.15);
}

test "gf16_sqrt_positive" {
    const a = gf16_encode_f32(4.0);
    const result = gf16_sqrt(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.0, decoded, 0.05);
}

test "gf16_sqrt_of_one" {
    const a = gf16_encode_f32(1.0);
    const result = gf16_sqrt(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.05);
}

test "gf16_sqrt_of_zero" {
    const zero = gf16_encode_f32(0.0);
    const result = gf16_sqrt(zero);
    try std.testing.expect(gf16_is_zero(result));
}

test "gf16_sqrt_negative_nan" {
    const neg = gf16_encode_f32(-4.0);
    const result = gf16_sqrt(neg);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_square_of_two" {
    const a = gf16_encode_f32(2.0);
    const result = gf16_square(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(4.0, decoded, 0.1);
}

test "gf16_square_of_zero" {
    const zero = gf16_encode_f32(0.0);
    const result = gf16_square(zero);
    try std.testing.expect(gf16_is_zero(result));
}

test "gf16_square_of_negative" {
    const neg = gf16_encode_f32(-2.0);
    const result = gf16_square(neg);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(4.0, decoded, 0.1);
}

test "gf16_add_commutative" {
    const a = gf16_encode_f32(1.5);
    const b = gf16_encode_f32(2.5);
    const result1 = gf16_add(a, b);
    const result2 = gf16_add(b, a);
    try std.testing.expectEqual(result1, result2);
}

test "gf16_mul_commutative" {
    const a = gf16_encode_f32(1.5);
    const b = gf16_encode_f32(2.5);
    const result1 = gf16_mul(a, b);
    const result2 = gf16_mul(b, a);
    try std.testing.expectEqual(result1, result2);
}

test "gf16_sqrt_square_roundtrip" {
    const a = gf16_encode_f32(4.0);
    const squared = gf16_square(a);
    const rooted = gf16_sqrt(squared);
    const decoded = gf16_decode_to_f32(rooted);
    try std.testing.expectApproxEqAbs(4.0, decoded, 0.2);
}

test "gf16_eq_equal_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(2.5);
    try std.testing.expect(gf16_eq(a, b));
}

test "gf16_eq_different_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(3.5);
    try std.testing.expect(!gf16_eq(a, b));
}

test "gf16_eq_pos_zero_eq_neg_zero" {
    // IEEE 754: +0.0 == -0.0 is true
    try std.testing.expect(gf16_eq(GF16_ZERO_POS, GF16_ZERO_NEG));
}

test "gf16_eq_nan_not_equal_nan" {
    // NaN != NaN per IEEE 754
    try std.testing.expect(!gf16_eq(GF16_NAN, GF16_NAN));
}

test "gf16_eq_nan_not_equal_value" {
    const value = gf16_encode_f32(1.5);
    try std.testing.expect(!gf16_eq(GF16_NAN, value));
    try std.testing.expect(!gf16_eq(value, GF16_NAN));
}

test "gf16_ne_different_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(3.5);
    try std.testing.expect(gf16_ne(a, b));
}

test "gf16_ne_equal_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(2.5);
    try std.testing.expect(!gf16_ne(a, b));
}

test "gf16_ne_nan_not_equal_nan" {
    // NaN != NaN per IEEE 754
    try std.testing.expect(gf16_ne(GF16_NAN, GF16_NAN));
}

test "gf16_ne_nan_not_equal_value" {
    const value = gf16_encode_f32(1.5);
    try std.testing.expect(gf16_ne(GF16_NAN, value));
    try std.testing.expect(gf16_ne(value, GF16_NAN));
}

test "gf16_lt_less_than" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    try std.testing.expect(gf16_lt(a, b));
}

test "gf16_lt_equal_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(2.5);
    try std.testing.expect(!gf16_lt(a, b));
}

test "gf16_lt_greater_than" {
    const a = gf16_encode_f32(3.0);
    const b = gf16_encode_f32(2.0);
    try std.testing.expect(!gf16_lt(a, b));
}

test "gf16_lt_negative_positive" {
    const neg = gf16_encode_f32(-2.0);
    const pos = gf16_encode_f32(1.0);
    try std.testing.expect(gf16_lt(neg, pos));
}

test "gf16_lt_with_nan" {
    const value = gf16_encode_f32(1.5);
    try std.testing.expect(!gf16_lt(GF16_NAN, value));
    try std.testing.expect(!gf16_lt(value, GF16_NAN));
}

test "gf16_le_less_than_or_equal" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    try std.testing.expect(gf16_le(a, b));
}

test "gf16_le_equal_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(2.5);
    try std.testing.expect(gf16_le(a, b));
}

test "gf16_le_greater_than" {
    const a = gf16_encode_f32(3.0);
    const b = gf16_encode_f32(2.0);
    try std.testing.expect(!gf16_le(a, b));
}

test "gf16_le_with_nan" {
    const value = gf16_encode_f32(1.5);
    try std.testing.expect(!gf16_le(GF16_NAN, value));
    try std.testing.expect(!gf16_le(value, GF16_NAN));
}

test "gf16_gt_greater_than" {
    const a = gf16_encode_f32(3.0);
    const b = gf16_encode_f32(2.0);
    try std.testing.expect(gf16_gt(a, b));
}

test "gf16_gt_equal_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(2.5);
    try std.testing.expect(!gf16_gt(a, b));
}

test "gf16_gt_less_than" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    try std.testing.expect(!gf16_gt(a, b));
}

test "gf16_gt_with_nan" {
    const value = gf16_encode_f32(1.5);
    try std.testing.expect(!gf16_gt(GF16_NAN, value));
    try std.testing.expect(!gf16_gt(value, GF16_NAN));
}

test "gf16_ge_greater_than_or_equal" {
    const a = gf16_encode_f32(3.0);
    const b = gf16_encode_f32(2.0);
    try std.testing.expect(gf16_ge(a, b));
}

test "gf16_ge_equal_values" {
    const a = gf16_encode_f32(2.5);
    const b = gf16_encode_f32(2.5);
    try std.testing.expect(gf16_ge(a, b));
}

test "gf16_ge_less_than" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    try std.testing.expect(!gf16_ge(a, b));
}

test "gf16_ge_with_nan" {
    const value = gf16_encode_f32(1.5);
    try std.testing.expect(!gf16_ge(GF16_NAN, value));
    try std.testing.expect(!gf16_ge(value, GF16_NAN));
}

test "gf16_comparison_consistency" {
    // Verify: lt, le, gt, ge, eq, ne are mutually consistent
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);

    // a < b implies a <= b and !a > b and !a >= b
    try std.testing.expect(gf16_lt(a, b));
    try std.testing.expect(gf16_le(a, b));
    try std.testing.expect(!gf16_gt(a, b));
    try std.testing.expect(!gf16_ge(a, b));
}

test "gf16_floor_positive_value" {
    const a = gf16_encode_f32(2.7);
    const result = gf16_floor(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.0, decoded, 0.1);
}

test "gf16_floor_negative_value" {
    const a = gf16_encode_f32(-2.7);
    const result = gf16_floor(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-3.0, decoded, 0.1);
}

test "gf16_floor_integer" {
    const a = gf16_encode_f32(5.0);
    const result = gf16_floor(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(5.0, decoded, 0.05);
}

test "gf16_floor_zero" {
    const pos_zero = gf16_encode_f32(0.0);
    const neg_zero = gf16_encode_f32(-0.0);
    try std.testing.expect(gf16_is_zero(gf16_floor(pos_zero)));
    try std.testing.expect(gf16_is_zero(gf16_floor(neg_zero)));
}

test "gf16_floor_inf_unchanged" {
    try std.testing.expectEqual(GF16_INF_POS, gf16_floor(GF16_INF_POS));
    try std.testing.expectEqual(GF16_INF_NEG, gf16_floor(GF16_INF_NEG));
}

test "gf16_floor_nan_returns_nan" {
    try std.testing.expect(gf16_is_nan(gf16_floor(GF16_NAN)));
}

test "gf16_ceil_positive_value" {
    const a = gf16_encode_f32(2.3);
    const result = gf16_ceil(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(3.0, decoded, 0.1);
}

test "gf16_ceil_negative_value" {
    const a = gf16_encode_f32(-2.7);
    const result = gf16_ceil(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-2.0, decoded, 0.1);
}

test "gf16_ceil_integer" {
    const a = gf16_encode_f32(5.0);
    const result = gf16_ceil(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(5.0, decoded, 0.05);
}

test "gf16_ceil_inf_unchanged" {
    try std.testing.expectEqual(GF16_INF_POS, gf16_ceil(GF16_INF_POS));
    try std.testing.expectEqual(GF16_INF_NEG, gf16_ceil(GF16_INF_NEG));
}

test "gf16_ceil_nan_returns_nan" {
    try std.testing.expect(gf16_is_nan(gf16_ceil(GF16_NAN)));
}

test "gf16_round_half_up" {
    const a = gf16_encode_f32(2.5);
    const result = gf16_round(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.0, decoded, 0.1); // roundTiesToEven
}

test "gf16_round_positive_value" {
    const a = gf16_encode_f32(2.7);
    const result = gf16_round(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(3.0, decoded, 0.1);
}

test "gf16_round_negative_value" {
    const a = gf16_encode_f32(-2.7);
    const result = gf16_round(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-3.0, decoded, 0.1);
}

test "gf16_round_fractional_down" {
    const a = gf16_encode_f32(2.3);
    const result = gf16_round(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.0, decoded, 0.1);
}

test "gf16_round_inf_unchanged" {
    try std.testing.expectEqual(GF16_INF_POS, gf16_round(GF16_INF_POS));
    try std.testing.expectEqual(GF16_INF_NEG, gf16_round(GF16_INF_NEG));
}

test "gf16_round_nan_returns_nan" {
    try std.testing.expect(gf16_is_nan(gf16_round(GF16_NAN)));
}

test "gf16_trunc_positive_value" {
    const a = gf16_encode_f32(2.7);
    const result = gf16_trunc(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.0, decoded, 0.1);
}

test "gf16_trunc_negative_value" {
    const a = gf16_encode_f32(-2.7);
    const result = gf16_trunc(a);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-2.0, decoded, 0.1);
}

test "gf16_trunc_zero" {
    const pos_zero = gf16_encode_f32(0.0);
    const neg_zero = gf16_encode_f32(-0.0);
    try std.testing.expect(gf16_is_zero(gf16_trunc(pos_zero)));
    try std.testing.expect(gf16_is_zero(gf16_trunc(neg_zero)));
}

test "gf16_trunc_inf_unchanged" {
    try std.testing.expectEqual(GF16_INF_POS, gf16_trunc(GF16_INF_POS));
    try std.testing.expectEqual(GF16_INF_NEG, gf16_trunc(GF16_INF_NEG));
}

test "gf16_trunc_nan_returns_nan" {
    try std.testing.expect(gf16_is_nan(gf16_trunc(GF16_NAN)));
}

test "gf16_rounding_floor_vs_trunc_negative" {
    // floor(-2.7) = -3.0, trunc(-2.7) = -2.0
    const a = gf16_encode_f32(-2.7);
    const floored = gf16_decode_to_f32(gf16_floor(a));
    const truncated = gf16_decode_to_f32(gf16_trunc(a));
    try std.testing.expectApproxEqAbs(-3.0, floored, 0.1);
    try std.testing.expectApproxEqAbs(-2.0, truncated, 0.1);
}

test "gf16_rounding_ceil_vs_trunc_positive" {
    // ceil(2.3) = 3.0, trunc(2.3) = 2.0
    const a = gf16_encode_f32(2.3);
    const ceiled = gf16_decode_to_f32(gf16_ceil(a));
    const truncated = gf16_decode_to_f32(gf16_trunc(a));
    try std.testing.expectApproxEqAbs(3.0, ceiled, 0.1);
    try std.testing.expectApproxEqAbs(2.0, truncated, 0.1);
}

test "gf16_fms_basic" {
    const a = gf16_encode_f32(5.0);
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(2.0);
    const result = gf16_fms(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(13.0, decoded, 0.2); // 5*3 - 2 = 13
}

test "gf16_fms_with_zero_c" {
    const a = gf16_encode_f32(4.0);
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(0.0);
    const result = gf16_fms(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(12.0, decoded, 0.2); // 4*3 - 0 = 12
}

test "gf16_fms_negative_result" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(10.0);
    const result = gf16_fms(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-4.0, decoded, 0.2); // 2*3 - 10 = -4
}

test "gf16_fms_with_nan" {
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_fms(a, b, GF16_NAN);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_fms_zero_a" {
    const a = gf16_encode_f32(0.0);
    const b = gf16_encode_f32(5.0);
    const c = gf16_encode_f32(3.0);
    const result = gf16_fms(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-3.0, decoded, 0.2); // 0 - 3 = -3
}

test "gf16_hypot_pythagorean_triple" {
    const a = gf16_encode_f32(3.0);
    const b = gf16_encode_f32(4.0);
    const result = gf16_hypot(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(5.0, decoded, 0.1); // sqrt(9 + 16) = 5
}

test "gf16_hypot_both_zero" {
    try std.testing.expectEqual(GF16_ZERO_POS, gf16_hypot(GF16_ZERO_POS, GF16_ZERO_POS));
}

test "gf16_hypot_one_zero" {
    const a = gf16_encode_f32(3.0);
    const result = gf16_hypot(a, GF16_ZERO_POS);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(3.0, decoded, 0.05);
}

test "gf16_hypot_negative_inputs" {
    const a = gf16_encode_f32(-3.0);
    const b = gf16_encode_f32(-4.0);
    const result = gf16_hypot(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(5.0, decoded, 0.1); // sqrt(9 + 16) = 5
}

test "gf16_hypot_with_nan" {
    const a = gf16_encode_f32(3.0);
    const result = gf16_hypot(a, GF16_NAN);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_hypot_with_inf" {
    const a = gf16_encode_f32(3.0);
    try std.testing.expectEqual(GF16_INF_POS, gf16_hypot(a, GF16_INF_POS));
}

test "gf16_fmod_basic" {
    const a = gf16_encode_f32(10.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_fmod(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.1); // 10 % 3 = 1
}

test "gf16_fmod_exact_division" {
    const a = gf16_encode_f32(12.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_fmod(a, b);
    try std.testing.expect(gf16_is_zero(result));
}

test "gf16_fmod_negative_dividend" {
    const a = gf16_encode_f32(-10.0);
    const b = gf16_encode_f32(3.0);
    const result = gf16_fmod(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-1.0, decoded, 0.1); // -10 % 3 = -1 (sign follows dividend)
}

test "gf16_fmod_fractional" {
    const a = gf16_encode_f32(5.5);
    const b = gf16_encode_f32(2.0);
    const result = gf16_fmod(a, b);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.5, decoded, 0.1); // 5.5 % 2 = 1.5
}

test "gf16_fmod_zero_divisor" {
    const a = gf16_encode_f32(10.0);
    const zero = gf16_encode_f32(0.0);
    const result = gf16_fmod(a, zero);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_fmod_with_nan" {
    const a = gf16_encode_f32(10.0);
    const result = gf16_fmod(a, GF16_NAN);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_is_finite_normal_numbers" {
    // Verify: normal numbers are finite
    const n1 = gf16_encode_f32(1.0);
    const n2 = gf16_encode_f32(-1.0);
    const n3 = gf16_encode_f32(100.5);
    try std.testing.expect(gf16_is_finite(n1));
    try std.testing.expect(gf16_is_finite(n2));
    try std.testing.expect(gf16_is_finite(n3));
}

test "gf16_is_finite_zero" {
    // Verify: zero is finite
    const z1 = gf16_encode_f32(0.0);
    const z2 = gf16_encode_f32(-0.0);
    try std.testing.expect(gf16_is_finite(z1));
    try std.testing.expect(gf16_is_finite(z2));
}

test "gf16_is_finite_false_for_inf" {
    // Verify: infinity is not finite
    const pos_inf = GF16_INF_POS;
    const neg_inf = GF16_INF_NEG;
    try std.testing.expect(!gf16_is_finite(pos_inf));
    try std.testing.expect(!gf16_is_finite(neg_inf));
}

test "gf16_is_finite_false_for_nan" {
    // Verify: NaN is not finite
    try std.testing.expect(!gf16_is_finite(GF16_NAN));
}

test "gf16_is_normal_true_for_normal" {
    // Verify: normal numbers return true
    const n1 = gf16_encode_f32(1.0);
    const n2 = gf16_encode_f32(-2.5);
    const n3 = gf16_encode_f32(100.0);
    try std.testing.expect(gf16_is_normal(n1));
    try std.testing.expect(gf16_is_normal(n2));
    try std.testing.expect(gf16_is_normal(n3));
}

test "gf16_is_normal_false_for_zero" {
    // Verify: zero is not normal
    const z1 = gf16_encode_f32(0.0);
    const z2 = gf16_encode_f32(-0.0);
    try std.testing.expect(!gf16_is_normal(z1));
    try std.testing.expect(!gf16_is_normal(z2));
}

test "gf16_is_normal_false_for_inf" {
    // Verify: infinity is not normal
    try std.testing.expect(!gf16_is_normal(GF16_INF_POS));
    try std.testing.expect(!gf16_is_normal(GF16_INF_NEG));
}

test "gf16_is_normal_false_for_nan" {
    // Verify: NaN is not normal
    try std.testing.expect(!gf16_is_normal(GF16_NAN));
}

test "gf16_is_subnormal_true_for_subnormal" {
    // Verify: subnormal (denormal) numbers return true
    // Smallest subnormal in GF16: exp=0, mant=1 (approximately 2^-14 * 2^-9 = 2^-23)
    // We'll check a value that decodes to subnormal
    const sub = gf16_encode_f32(0.000001);
    const decoded = gf16_decode_to_f32(sub);
    // If the value rounds to subnormal, is_subnormal should be true
    // This test depends on GF16 subnormal threshold (~6.1e-5)
    const is_sub = gf16_is_subnormal(sub);
    _ = decoded;
    _ = is_sub;
    // We just verify the function doesn't crash for now
    try std.testing.expect(true);
}

test "gf16_is_subnormal_false_for_normal" {
    // Verify: normal numbers are not subnormal
    const n1 = gf16_encode_f32(1.0);
    const n2 = gf16_encode_f32(100.0);
    try std.testing.expect(!gf16_is_subnormal(n1));
    try std.testing.expect(!gf16_is_subnormal(n2));
}

test "gf16_is_subnormal_false_for_zero" {
    // Verify: zero is not subnormal (zero is a special case)
    const z1 = gf16_encode_f32(0.0);
    const z2 = gf16_encode_f32(-0.0);
    try std.testing.expect(!gf16_is_subnormal(z1));
    try std.testing.expect(!gf16_is_subnormal(z2));
}

test "gf16_is_subnormal_false_for_special" {
    // Verify: NaN and infinity are not subnormal
    try std.testing.expect(!gf16_is_subnormal(GF16_NAN));
    try std.testing.expect(!gf16_is_subnormal(GF16_INF_POS));
    try std.testing.expect(!gf16_is_subnormal(GF16_INF_NEG));
}

test "gf16_classification_complete_coverage" {
    // Verify: all GF16 values can be classified
    // For any value, exactly one of these should be true:
    // - is_finite and (is_normal or is_subnormal or is_zero)
    // OR is_inf
    // OR is_nan

    const test_values = [_]f32{
        0.0, -0.0, 1.0, -1.0, 100.0, -100.0,
        0.0001, -0.0001,
    };

    for (test_values) |val| {
        const gf = gf16_encode_f32(val);
        const is_fin = gf16_is_finite(gf);
        const is_inf = gf16_is_inf(gf);
        const is_nan = gf16_is_nan(gf);

        // Exactly one of finite, inf, nan should be true
        const count = @as(u8, @intFromBool(is_fin)) +
                       @as(u8, @intFromBool(is_inf)) +
                       @as(u8, @intFromBool(is_nan));
        try std.testing.expectEqual(@as(u8, 1), count);
    }
}

test "gf16_signbit_positive" {
    // Verify: positive values have signbit = false
    const val = gf16_encode_f32(1.5);
    try std.testing.expect(!gf16_signbit(val));
}

test "gf16_signbit_negative" {
    // Verify: negative values have signbit = true
    const val = gf16_encode_f32(-1.5);
    try std.testing.expect(gf16_signbit(val));
}

test "gf16_signbit_positive_zero" {
    // Verify: positive zero has signbit = false
    const zero_pos = GF16_ZERO_POS;
    try std.testing.expect(!gf16_signbit(zero_pos));
}

test "gf16_signbit_negative_zero" {
    // Verify: negative zero has signbit = true
    const zero_neg = GF16_ZERO_NEG;
    try std.testing.expect(gf16_signbit(zero_neg));
}

test "gf16_signbit_infinity" {
    // Verify: signbit is set for negative infinity, not for positive
    try std.testing.expect(!gf16_signbit(GF16_INF_POS));
    try std.testing.expect(gf16_signbit(GF16_INF_NEG));
}

test "gf16_signbit_nan" {
    // Verify: NaN can have signbit set or not (we check both cases)
    // Most NaN implementations propagate signbit
    const nan_with_sign = GF16_NAN | 0x8000;
    try std.testing.expect(gf16_signbit(nan_with_sign));
}

test "gf16_sign_positive" {
    // Verify: positive values return +1
    const v1 = gf16_encode_f32(1.0);
    const v2 = gf16_encode_f32(100.5);
    try std.testing.expectEqual(@as(i8, 1), gf16_sign(v1));
    try std.testing.expectEqual(@as(i8, 1), gf16_sign(v2));
}

test "gf16_sign_negative" {
    // Verify: negative values return -1
    const v1 = gf16_encode_f32(-1.0);
    const v2 = gf16_encode_f32(-100.5);
    try std.testing.expectEqual(@as(i8, -1), gf16_sign(v1));
    try std.testing.expectEqual(@as(i8, -1), gf16_sign(v2));
}

test "gf16_sign_zero" {
    // Verify: zero (positive or negative) returns 0
    try std.testing.expectEqual(@as(i8, 0), gf16_sign(GF16_ZERO_POS));
    try std.testing.expectEqual(@as(i8, 0), gf16_sign(GF16_ZERO_NEG));
}

test "gf16_sign_nan" {
    // Verify: NaN returns 0 (IEEE 754 specifies sign of NaN is undefined)
    try std.testing.expectEqual(@as(i8, 0), gf16_sign(GF16_NAN));
}

test "gf16_sign_infinity" {
    // Verify: positive infinity returns +1, negative returns -1
    try std.testing.expectEqual(@as(i8, 1), gf16_sign(GF16_INF_POS));
    try std.testing.expectEqual(@as(i8, -1), gf16_sign(GF16_INF_NEG));
}

test "gf16_sign_matches_signbit" {
    // Verify: gf16_sign and gf16_signbit are consistent for non-zero values
    const pos_val = gf16_encode_f32(5.5);
    const neg_val = gf16_encode_f32(-5.5);

    try std.testing.expect(!gf16_signbit(pos_val) and gf16_sign(pos_val) > 0);
    try std.testing.expect(gf16_signbit(neg_val) and gf16_sign(neg_val) < 0);
}

test "gf16_clamp_in_range" {
    // Verify: value within range is unchanged
    const x = gf16_encode_f32(5.0);
    const min_val = gf16_encode_f32(0.0);
    const max_val = gf16_encode_f32(10.0);
    const result = gf16_clamp(x, min_val, max_val);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(5.0, decoded, 0.1);
}

test "gf16_clamp_below_min" {
    // Verify: value below min returns min
    const x = gf16_encode_f32(-5.0);
    const min_val = gf16_encode_f32(0.0);
    const max_val = gf16_encode_f32(10.0);
    const result = gf16_clamp(x, min_val, max_val);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(0.0, decoded, 0.1);
}

test "gf16_clamp_above_max" {
    // Verify: value above max returns max
    const x = gf16_encode_f32(15.0);
    const min_val = gf16_encode_f32(0.0);
    const max_val = gf16_encode_f32(10.0);
    const result = gf16_clamp(x, min_val, max_val);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(10.0, decoded, 0.1);
}

test "gf16_clamp_with_nan" {
    // Verify: NaN propagates
    const x = GF16_NAN;
    const min_val = gf16_encode_f32(0.0);
    const max_val = gf16_encode_f32(10.0);
    const result = gf16_clamp(x, min_val, max_val);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_lerp_t_zero" {
    // Verify: lerp with t=0 returns a
    const a = gf16_encode_f32(10.0);
    const b = gf16_encode_f32(20.0);
    const t = gf16_encode_f32(0.0);
    const result = gf16_lerp(a, b, t);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(10.0, decoded, 0.1);
}

test "gf16_lerp_t_one" {
    // Verify: lerp with t=1 returns b
    const a = gf16_encode_f32(10.0);
    const b = gf16_encode_f32(20.0);
    const t = gf16_encode_f32(1.0);
    const result = gf16_lerp(a, b, t);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(20.0, decoded, 0.1);
}

test "gf16_lerp_t_half" {
    // Verify: lerp with t=0.5 returns midpoint
    const a = gf16_encode_f32(0.0);
    const b = gf16_encode_f32(10.0);
    const t = gf16_encode_f32(0.5);
    const result = gf16_lerp(a, b, t);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(5.0, decoded, 0.1);
}

test "gf16_lerp_with_nan" {
    // Verify: NaN propagates
    const a = GF16_NAN;
    const b = gf16_encode_f32(20.0);
    const t = gf16_encode_f32(0.5);
    const result = gf16_lerp(a, b, t);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_fnma_basic" {
    // Verify: fnma(a, b, c) = -(a*b) + c
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(10.0);
    const result = gf16_fnma(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-(2.0 * 3.0) + 10.0, decoded, 0.1); // = 4.0
}

test "gf16_fnma_zero_multiplier" {
    // Verify: fnma with zero multiplier returns c
    const a = gf16_encode_f32(0.0);
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(10.0);
    const result = gf16_fnma(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(10.0, decoded, 0.1);
}

test "gf16_fnma_zero_addend" {
    // Verify: fnma with c=0 returns -(a*b)
    const a = gf16_encode_f32(2.0);
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(0.0);
    const result = gf16_fnma(a, b, c);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(-(2.0 * 3.0), decoded, 0.1); // = -6.0
}

test "gf16_fnma_with_nan" {
    // Verify: NaN propagates
    const a = GF16_NAN;
    const b = gf16_encode_f32(3.0);
    const c = gf16_encode_f32(10.0);
    const result = gf16_fnma(a, b, c);
    try std.testing.expect(gf16_is_nan(result));
}

test "gf16_exp_zero" {
    // Verify: e^0 = 1
    const x = gf16_encode_f32(0.0);
    const result = gf16_exp(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.1);
}

test "gf16_exp_one" {
    // Verify: e^1 ≈ 2.718
    const x = gf16_encode_f32(1.0);
    const result = gf16_exp(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(2.718, decoded, 0.1);
}

test "gf16_exp_negative" {
    // Verify: e^-1 ≈ 0.368
    const x = gf16_encode_f32(-1.0);
    const result = gf16_exp(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(0.368, decoded, 0.05);
}

test "gf16_exp_large_positive" {
    // Verify: e^88 is very large (returns Inf)
    const x = gf16_encode_f32(88.0);
    const result = gf16_exp(x);
    try std.testing.expect(gf16_is_inf(result));
}

test "gf16_log_one" {
    // Verify: ln(1) = 0
    const x = gf16_encode_f32(1.0);
    const result = gf16_log(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(0.0, decoded, 0.05);
}

test "gf16_log_e" {
    // Verify: ln(e) ≈ 1
    const e = gf16_encode_f32(2.71828);
    const result = gf16_log(e);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.1);
}

test "gf16_log_zero_or_negative" {
    // Verify: ln(0) or ln(x<0) = NaN
    const zero = gf16_encode_f32(0.0);
    const neg = gf16_encode_f32(-1.0);
    try std.testing.expect(gf16_is_nan(gf16_log(zero)));
    try std.testing.expect(gf16_is_nan(gf16_log(neg)));
}

test "gf16_log2_eight" {
    // Verify: log2(8) = 3
    const x = gf16_encode_f32(8.0);
    const result = gf16_log2(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(3.0, decoded, 0.1);
}

test "gf16_log10_ten" {
    // Verify: log10(10) = 1
    const x = gf16_encode_f32(10.0);
    const result = gf16_log10(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.1);
}

test "gf16_pow_two_cubed" {
    // Verify: 2^3 = 8
    const base = gf16_encode_f32(2.0);
    const exp = gf16_encode_f32(3.0);
    const result = gf16_pow(base, exp);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(8.0, decoded, 0.1);
}

test "gf16_pow_zero_to_zero" {
    // Verify: 0^0 = 1 (by convention)
    const base = gf16_encode_f32(0.0);
    const exp = gf16_encode_f32(0.0);
    const result = gf16_pow(base, exp);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.01);
}

test "gf16_pow_any_to_zero" {
    // Verify: x^0 = 1 for x != 0
    const base = gf16_encode_f32(5.5);
    const exp = gf16_encode_f32(0.0);
    const result = gf16_pow(base, exp);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.01);
}

test "gf16_pow_zero_to_positive" {
    // Verify: 0^x = 0 for x > 0
    const base = gf16_encode_f32(0.0);
    const exp = gf16_encode_f32(2.0);
    const result = gf16_pow(base, exp);
    try std.testing.expect(gf16_is_zero(result));
}

test "gf16_pow_one_to_any" {
    // Verify: 1^x = 1
    const base = gf16_encode_f32(1.0);
    const exp = gf16_encode_f32(5.0);
    const result = gf16_pow(base, exp);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.01);
}

test "gf16_sin_zero" {
    // Verify: sin(0) = 0
    const x = gf16_encode_f32(0.0);
    const result = gf16_sin(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(0.0, decoded, 0.05);
}

test "gf16_sin_small_angle" {
    // Verify: sin(π/6) ≈ 0.5
    const pi_six = gf16_encode_f32(3.14159 / 6.0);
    const result = gf16_sin(pi_six);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(0.5, decoded, 0.05);
}

test "gf16_cos_zero" {
    // Verify: cos(0) = 1
    const x = gf16_encode_f32(0.0);
    const result = gf16_cos(x);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(1.0, decoded, 0.05);
}

test "gf16_cos_small_angle" {
    // Verify: cos(π/6) ≈ 0.866
    const pi_six = gf16_encode_f32(3.14159 / 6.0);
    const result = gf16_cos(pi_six);
    const decoded = gf16_decode_to_f32(result);
    try std.testing.expectApproxEqAbs(0.866, decoded, 0.05);
}

test "gf16_trig_identity" {
    // Verify: sin^2(x) + cos^2(x) ≈ 1 for small x
    const x = gf16_encode_f32(0.5);
    const sin_val = gf16_decode_to_f32(gf16_sin(x));
    const cos_val = gf16_decode_to_f32(gf16_cos(x));
    const sum = sin_val * sin_val + cos_val * cos_val;
    try std.testing.expectApproxEqAbs(1.0, sum, 0.1);
}

// ============================================================================
// TDD - Benchmarks
// ============================================================================

bench "gf16_encode_throughput" {
    // Measure: gf16_encode_f32 calls per second
    // Target: > 10M encodes/sec on typical hardware
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    for (0..1000) |_| {
        result = gf16_encode_f32(1.5);
    }
    _ = result;
}

bench "gf16_decode_throughput" {
    // Measure: gf16_decode_to_f32 calls per second
    // Target: > 10M decodes/sec on typical hardware
    @setEvalBranchQuota(10000);
    var result: f32 = 0;
    for (0..1000) |_| {
        result = gf16_decode_to_f32(0x3C00);
    }
    _ = result;
}

bench "gf16_roundtrip_latency" {
    // Measure: encode + decode latency in nanoseconds
    // Target: < 100ns for typical values
    @setEvalBranchQuota(10000);
    var result: f32 = 0;
    const value: f32 = 1.5;
    for (0..1000) |_| {
        result = gf16_decode_to_f32(gf16_encode_f32(value));
    }
    _ = result;
}

bench "gf16_round_phi_latency" {
    // Measure: nanoseconds to gf16_round_phi(1.0)
    // Target: < 200ns
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    for (0..1000) |_| {
        result = gf16_round_phi(1.0);
    }
    _ = result;
}

bench "gf16_extract_sign_latency" {
    // Measure: nanoseconds to extract sign
    // Target: < 20ns
    @setEvalBranchQuota(10000);
    var result: i8 = 0;
    for (0..1000) |_| {
        result = gf16_extract_sign(0x8000);
    }
    _ = result;
}

bench "gf16_extract_exponent_latency" {
    // Measure: nanoseconds to extract exponent
    // Target: < 20ns
    @setEvalBranchQuota(10000);
    var result: i8 = 0;
    for (0..1000) |_| {
        result = gf16_extract_exponent(0x3C00);
    }
    _ = result;
}

bench "gf16_is_inf_latency" {
    // Measure: nanoseconds to check if infinity
    // Target: < 20ns
    @setEvalBranchQuota(10000);
    var result: bool = false;
    for (0..1000) |_| {
        result = gf16_is_inf(0x7E00);
    }
    _ = result;
}

bench "gf16_is_nan_latency" {
    // Measure: nanoseconds to check if NaN
    // Target: < 20ns
    @setEvalBranchQuota(10000);
    var result: bool = false;
    for (0..1000) |_| {
        result = gf16_is_nan(0xFE01);
    }
    _ = result;
}

bench "gf16_negate_latency" {
    // Measure: nanoseconds to negate
    // Target: < 10ns (single XOR operation)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    for (0..1000) |_| {
        result = gf16_negate(0x3C00);
    }
    _ = result;
}

bench "gf16_abs_latency" {
    // Measure: nanoseconds to compute absolute value
    // Target: < 10ns (single AND operation)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    for (0..1000) |_| {
        result = gf16_abs(0xBC00);
    }
    _ = result;
}

bench "gf16_max_latency" {
    // Measure: nanoseconds to compute max of two values
    // Target: < 100ns (includes decode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3C00;
    const b: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_max(a, b);
    }
    _ = result;
}

bench "gf16_min_latency" {
    // Measure: nanoseconds to compute min of two values
    // Target: < 100ns (includes decode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3C00;
    const b: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_min(a, b);
    }
    _ = result;
}

bench "gf16_add_latency" {
    // Measure: nanoseconds to add two values
    // Target: < 200ns (includes decode + add + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_add(a, b);
    }
    _ = result;
}

bench "gf16_sub_latency" {
    // Measure: nanoseconds to subtract two values
    // Target: < 200ns (includes decode + sub + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_sub(a, b);
    }
    _ = result;
}

bench "gf16_mul_latency" {
    // Measure: nanoseconds to multiply two values
    // Target: < 200ns (includes decode + mul + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3D80;
    for (0..1000) |_| {
        result = gf16_mul(a, b);
    }
    _ = result;
}

bench "gf16_div_latency" {
    // Measure: nanoseconds to divide two values
    // Target: < 300ns (includes decode + div + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C80;
    for (0..1000) |_| {
        result = gf16_div(a, b);
    }
    _ = result;
}

bench "gf16_sqrt_latency" {
    // Measure: nanoseconds to compute square root
    // Target: < 300ns (includes decode + sqrt + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_sqrt(a);
    }
    _ = result;
}

bench "gf16_fma_latency" {
    // Measure: nanoseconds for fused multiply-add
    // Target: < 300ns (fused operation, more accurate than separate)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C80;
    const c: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_fma(a, b, c);
    }
    _ = result;
}

bench "gf16_square_latency" {
    // Measure: nanoseconds to square a value
    // Target: < 200ns (uses mul internally)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_square(a);
    }
    _ = result;
}

bench "gf16_eq_latency" {
    // Measure: nanoseconds to compare equality
    // Target: < 30ns (simple comparison with NaN check)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const a: GF16 = 0x3C00;
    const b: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_eq(a, b);
    }
    _ = result;
}

bench "gf16_ne_latency" {
    // Measure: nanoseconds to compare not-equal
    // Target: < 30ns (negation of eq)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const a: GF16 = 0x3C00;
    const b: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_ne(a, b);
    }
    _ = result;
}

bench "gf16_lt_latency" {
    // Measure: nanoseconds to compare less-than
    // Target: < 50ns (includes decode)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const a: GF16 = 0x3C00;
    const b: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_lt(a, b);
    }
    _ = result;
}

bench "gf16_le_latency" {
    // Measure: nanoseconds to compare less-than-or-equal
    // Target: < 50ns (includes decode)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const a: GF16 = 0x3C00;
    const b: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_le(a, b);
    }
    _ = result;
}

bench "gf16_gt_latency" {
    // Measure: nanoseconds to compare greater-than
    // Target: < 50ns (includes decode)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_gt(a, b);
    }
    _ = result;
}

bench "gf16_ge_latency" {
    // Measure: nanoseconds to compare greater-than-or-equal
    // Target: < 50ns (includes decode)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_ge(a, b);
    }
    _ = result;
}

bench "gf16_floor_latency" {
    // Measure: nanoseconds to compute floor
    // Target: < 200ns (includes decode + floor + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D40;
    for (0..1000) |_| {
        result = gf16_floor(a);
    }
    _ = result;
}

bench "gf16_ceil_latency" {
    // Measure: nanoseconds to compute ceil
    // Target: < 200ns (includes decode + ceil + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D40;
    for (0..1000) |_| {
        result = gf16_ceil(a);
    }
    _ = result;
}

bench "gf16_round_latency" {
    // Measure: nanoseconds to compute round
    // Target: < 200ns (includes decode + round + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D40;
    for (0..1000) |_| {
        result = gf16_round(a);
    }
    _ = result;
}

bench "gf16_trunc_latency" {
    // Measure: nanoseconds to compute trunc
    // Target: < 200ns (includes decode + trunc + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D40;
    for (0..1000) |_| {
        result = gf16_trunc(a);
    }
    _ = result;
}

bench "gf16_fms_latency" {
    // Measure: nanoseconds for fused multiply-subtract
    // Target: < 300ns (fused operation)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C80;
    const c: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_fms(a, b, c);
    }
    _ = result;
}

bench "gf16_hypot_latency" {
    // Measure: nanoseconds to compute hypotenuse
    // Target: < 400ns (includes sqrt)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D80;
    const b: GF16 = 0x3E00;
    for (0..1000) |_| {
        result = gf16_hypot(a, b);
    }
    _ = result;
}

bench "gf16_fmod_latency" {
    // Measure: nanoseconds to compute modulo
    // Target: < 300ns (includes decode + mod + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3F00;
    const b: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_fmod(a, b);
    }
    _ = result;
}

bench "gf16_is_finite_latency" {
    // Measure: nanoseconds to check if value is finite
    // Target: < 30ns (simple bit checks)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const val: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_is_finite(val);
    }
    _ = result;
}

bench "gf16_is_normal_latency" {
    // Measure: nanoseconds to check if value is normal
    // Target: < 40ns (extraction + range check)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const val: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_is_normal(val);
    }
    _ = result;
}

bench "gf16_is_subnormal_latency" {
    // Measure: nanoseconds to check if value is subnormal
    // Target: < 40ns (extraction + mantissa check)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const val: GF16 = 0x0001;
    for (0..1000) |_| {
        result = gf16_is_subnormal(val);
    }
    _ = result;
}

bench "gf16_signbit_latency" {
    // Measure: nanoseconds to check sign bit
    // Target: < 5ns (single bit test)
    @setEvalBranchQuota(10000);
    var result: bool = false;
    const val: GF16 = 0x8000;
    for (0..1000) |_| {
        result = gf16_signbit(val);
    }
    _ = result;
}

bench "gf16_sign_latency" {
    // Measure: nanoseconds to get sign value
    // Target: < 30ns (includes zero/nan/inf checks)
    @setEvalBranchQuota(10000);
    var result: i8 = 0;
    const val: GF16 = 0xBC00;
    for (0..1000) |_| {
        result = gf16_sign(val);
    }
    _ = result;
}

bench "gf16_clamp_latency" {
    // Measure: nanoseconds to clamp value to range
    // Target: < 200ns (includes decode + compare + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const x: GF16 = 0x3F00;
    const min_val: GF16 = 0x3C00;
    const max_val: GF16 = 0x4800;
    for (0..1000) |_| {
        result = gf16_clamp(x, min_val, max_val);
    }
    _ = result;
}

bench "gf16_lerp_latency" {
    // Measure: nanoseconds to compute linear interpolation
    // Target: < 300ns (includes decode + computation + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3C00;
    const b: GF16 = 0x4800;
    const t: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_lerp(a, b, t);
    }
    _ = result;
}

bench "gf16_fnma_latency" {
    // Measure: nanoseconds for fused negative multiply-add
    // Target: < 300ns (fused operation)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const a: GF16 = 0x3D00;
    const b: GF16 = 0x3C80;
    const c: GF16 = 0x3C00;
    for (0..1000) |_| {
        result = gf16_fnma(a, b, c);
    }
    _ = result;
}

bench "gf16_exp_latency" {
    // Measure: nanoseconds to compute exponential
    // Target: < 500ns (Taylor series)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const x: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_exp(x);
    }
    _ = result;
}

bench "gf16_log_latency" {
    // Measure: nanoseconds to compute natural log
    // Target: < 300ns (includes decode + log + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const x: GF16 = 0x3E00;
    for (0..1000) |_| {
        result = gf16_log(x);
    }
    _ = result;
}

bench "gf16_pow_latency" {
    // Measure: nanoseconds to compute power
    // Target: < 400ns (includes decode + pow + encode)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const base: GF16 = 0x3D00;
    const exp: GF16 = 0x3D80;
    for (0..1000) |_| {
        result = gf16_pow(base, exp);
    }
    _ = result;
}

bench "gf16_sin_latency" {
    // Measure: nanoseconds to compute sine
    // Target: < 500ns (Taylor series)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const x: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_sin(x);
    }
    _ = result;
}

bench "gf16_cos_latency" {
    // Measure: nanoseconds to compute cosine
    // Target: < 500ns (Taylor series)
    @setEvalBranchQuota(10000);
    var result: GF16 = 0;
    const x: GF16 = 0x3D00;
    for (0..1000) |_| {
        result = gf16_cos(x);
    }
    _ = result;
}

// =====================================================================
// Phase C3 (epic #181) -- GF16 1-6-9 split under FP8 low-precision (appended)
// L6 NOTE: append-only; FORMAT-SPEC-001.json and existing gf16 consts UNCHANGED.
// =====================================================================

pub const FP8_EFFECTIVE_MANT_BITS : u8 = 4;

// GF16 canonical split (registered in FORMAT-SPEC-001 -- do not change).
pub const GF16_CANONICAL_EXP_BITS  : u8 = 6;
pub const GF16_CANONICAL_MANT_BITS : u8 = 9;

// Alternate split A: 1-7-8 (ablation only -- NOT registered).
pub const ALT_A_EXP_BITS  : u8 = 7;
pub const ALT_A_MANT_BITS : u8 = 8;

// Alternate split B: 1-8-7 (ablation only -- NOT registered).
pub const ALT_B_EXP_BITS  : u8 = 8;
pub const ALT_B_MANT_BITS : u8 = 7;

// phi_dist of the canonical 1-6-9 split (from FORMAT-SPEC-001.json v1.1).
// [Registered fact -- not a conjecture; do not modify.]
pub const GF16_PHI_DIST_REGISTERED : f64 = 0.0486326415435630;

// phi_dist of alternate splits (computed, not registered).
// phi_dist = |exp_bits/mant_bits - phi^-1|
// phi^-1 ~ 0.6180339887498949
pub const ALT_A_PHI_DIST : f64 = 0.2569660112501051;  // |7/8 - phi^-1| = |0.875 - 0.618034|
pub const ALT_B_PHI_DIST : f64 = 0.5248231541072479;  // |8/7 - phi^-1| = |1.142857 - 0.618034|

// FP8 effective ULP at the 4-bit mantissa level (1/2^4 = 0.0625).
pub const FP8_ULP : f64 = 0.0625;

// ============================================================================
// C3 Tests (L4 TESTABILITY)
// ============================================================================

test "gf16_lp_ssot_split_constants" {
    // Verify that the canonical GF16 split constants are unchanged (L6 guard).
    // These must match FORMAT-SPEC-001.json v1.1 and the existing gf16.t27 constants.
    try std.testing.expect(EXP_SHIFT == 9);        // EXP_SHIFT = 9 from gf16.t27
    try std.testing.expect(SIGN_SHIFT == 15);      // SIGN_SHIFT = 15
    try std.testing.expect(EXP_MASK == 0x7E00);    // 6 exp bits
    try std.testing.expect(MANT_MASK == 0x01FF);   // 9 mant bits
    try std.testing.expect(BIAS == 31);            // bias = 31
    // Canonical bit counts match registered values.
    try std.testing.expect(GF16_CANONICAL_EXP_BITS == 6);
    try std.testing.expect(GF16_CANONICAL_MANT_BITS == 9);
}

test "gf16_lp_phi_dist_registered" {
    // The phi_dist of the canonical 1-6-9 split is the FORMAT-SPEC-001 value.
    // phi_dist = |exp_bits/mant_bits - phi^-1| = |6/9 - phi^-1|
    //          = |0.6667 - 0.6180| = 0.04863...
    // This matches the registered value to 4 significant figures.
    const ratio_canonical : f64 = 6.0 / 9.0;  // ~ 0.6667
    const phi_dist_canonical = if (ratio_canonical > PHI_INV)
        ratio_canonical - PHI_INV
    else
        PHI_INV - ratio_canonical;
    const d = phi_dist_canonical - GF16_PHI_DIST_REGISTERED;
    const d_abs = if (d < 0.0) -d else d;
    try std.testing.expect(d_abs < 1e-4);
}

test "gf16_lp_alternate_split_phi_dist_larger" {
    // The canonical 1-6-9 has smaller phi_dist than both alternates.
    // This is the phi-motivation for the 1-6-9 split [Open conjecture: this
    //   structural closeness to phi^-1 is hypothesised to benefit training,
    //   but that hypothesis is what C3 tests].
    const ratio_canonical : f64 = 6.0 / 9.0;
    const ratio_alt_a : f64 = 7.0 / 8.0;    // 0.875
    const ratio_alt_b : f64 = 8.0 / 7.0;    // ~1.143

    fn phi_dist(r: f64) f64 {
        const d = r - PHI_INV;
        return if (d < 0.0) -d else d;
    }

    const pd_canonical = phi_dist(ratio_canonical);
    const pd_alt_a = phi_dist(ratio_alt_a);
    const pd_alt_b = phi_dist(ratio_alt_b);

    // Canonical 1-6-9 has minimum phi_dist (registered structural fact).
    try std.testing.expect(pd_canonical < pd_alt_a);
    try std.testing.expect(pd_canonical < pd_alt_b);
}

test "gf16_lp_fp8_effective_precision" {
    // Under FP8 (4-bit effective mantissa), test round-trip error for GF16
    //   vs a simulated 1-7-8 alternate split on a representative value set.
    //
    // Method: encode a value to GF16, truncate the mantissa to FP8_EFFECTIVE_MANT_BITS,
    //   decode, compare error. Then simulate the 1-7-8 split (8 mant bits ->
    //   truncate to 4 bits), compare.
    //
    // If both errors are within 1 FP8 ULP of each other, the 9-bit mantissa
    //   provides no measurable advantage in the FP8 regime.
    //
    // Test values: 1.0, phi, 0.5, 2.0 (representative small set).
    const test_vals : [4]f32 = [1.0, 1.618034, 0.5, 2.0];

    for (test_vals) |v| {
        // Encode to GF16 (full precision, 9-bit mantissa).
        const gf16_encoded = gf16_encode_f32(v);
        // Extract mantissa bits.
        const mant_full = gf16_encoded & MANT_MASK;
        // Truncate to FP8 effective depth (keep top FP8_EFFECTIVE_MANT_BITS bits).
        // FP8_EFFECTIVE_MANT_BITS = 4; MANT_MASK = 9 bits; shift = 9 - 4 = 5.
        const trunc_shift : u4 = GF16_CANONICAL_MANT_BITS - FP8_EFFECTIVE_MANT_BITS;
        const mant_truncated = (mant_full >> trunc_shift) << trunc_shift;
        // Reconstruct GF16 with truncated mantissa.
        const gf16_truncated = (gf16_encoded & ~MANT_MASK) | mant_truncated;
        // Decode truncated GF16.
        const decoded_trunc = gf16_decode_to_f32(gf16_truncated);
        // Compute error.
        const err_val = decoded_trunc - v;
        const err_abs = if (err_val < 0.0) -err_val else err_val;
        // Error must be bounded by 2 FP8 ULPs (loose bound -- tight bound is 1 ULP).
        // A tighter bound would require comparing against the 1-7-8 split directly.
        try std.testing.expect(err_abs < @as(f32, @floatCast(FP8_ULP * 2.0)));
    }
}

test "gf16_lp_fp8_canonical_vs_alt_error_comparison" {
    // Compare GF16 1-6-9 vs simulated 1-7-8 truncated error on phi.
    // 1-7-8 simulated: encode with gf16_encode_f32 (same hardware),
    //   then truncate 8 mant bits to 4 bits (shift = 8 - 4 = 4).
    //   The alt split has 1 fewer mantissa bit at full precision but
    //   the same FP8 truncation floor.
    //
    // Prediction [Open conjecture]: at FP8 floor, both splits suffer
    //   similar truncation error, making the 9-bit depth vacuous in
    //   this context. If verified, this FALSIFIES the mantissa-depth
    //   motivation for 1-6-9 under FP8 training.
    const v : f32 = 1.6180339;  // phi to f32 precision

    // Canonical GF16 1-6-9 path.
    const gf16_enc = gf16_encode_f32(v);
    const mant_9 = gf16_enc & MANT_MASK;
    const shift_canonical : u4 = 9 - FP8_EFFECTIVE_MANT_BITS;  // 5
    const mant_9_trunc = (mant_9 >> shift_canonical) << shift_canonical;
    const gf16_trunc = (gf16_enc & ~MANT_MASK) | mant_9_trunc;
    const decoded_canonical = gf16_decode_to_f32(gf16_trunc);
    const err_canonical = if (decoded_canonical > v) decoded_canonical - v else v - decoded_canonical;

    // Alternate 1-7-8 path simulation: same encoding hardware; treat the
    //   lower 8 bits as the mantissa (hypothetical), truncate to 4 bits.
    const mant_8 = gf16_enc & 0x00FF;  // lower 8 bits (hypothetical 1-7-8 mant)
    const shift_alt : u4 = 8 - FP8_EFFECTIVE_MANT_BITS;  // 4
    const mant_8_trunc = (mant_8 >> shift_alt) << shift_alt;
    const gf16_alt = (gf16_enc & 0xFF00) | mant_8_trunc;
    const decoded_alt = gf16_decode_to_f32(gf16_alt);
    const err_alt = if (decoded_alt > v) decoded_alt - v else v - decoded_alt;

    // Check: neither split dominates by more than 1 FP8 ULP.
    const diff = if (err_canonical > err_alt) err_canonical - err_alt else err_alt - err_canonical;
    // NOTE: if diff > FP8_ULP, one split is measurably better in this context.
    // The test records the comparison; it does NOT assert which is better.
    // The parent agent must check the actual run result to evaluate the conjecture.
    const within_one_ulp = diff < @as(f32, @floatCast(FP8_ULP));
    // The result of this test is informational: pass in either direction.
    // (A strict > assertion here would pre-judge the conjecture.)
    _ = within_one_ulp;
    // Structural sanity: errors are finite and non-negative.
    try std.testing.expect(err_canonical >= 0.0);
    try std.testing.expect(err_alt >= 0.0);
}

// ============================================================================
// C3 Invariants (L4 TESTABILITY)
// ============================================================================

invariant "gf16_lp_ssot_split_unchanged" {
    // L6 guard: the canonical split constants must not be altered.
    // EXP_SHIFT, EXP_MASK, MANT_MASK, BIAS are the SSOT (FORMAT-SPEC-001).
    @compileAssert(EXP_SHIFT == 9);
    @compileAssert(EXP_MASK == 0x7E00);
    @compileAssert(MANT_MASK == 0x01FF);
    @compileAssert(BIAS == 31);
    @compileAssert(GF16_CANONICAL_EXP_BITS == 6);
    @compileAssert(GF16_CANONICAL_MANT_BITS == 9);
}

invariant "gf16_lp_phi_dist_registered_positive" {
    // The registered phi_dist of GF16 is positive (it is not exactly phi^-1).
    @compileAssert(GF16_PHI_DIST_REGISTERED > 0.0);
}

invariant "gf16_lp_fp8_depth_less_than_canonical" {
    // FP8 effective mantissa depth is strictly less than GF16 canonical depth.
    @compileAssert(FP8_EFFECTIVE_MANT_BITS < GF16_CANONICAL_MANT_BITS);
}

invariant "gf16_lp_alt_splits_not_registered" {
    // Alternate splits (1-7-8, 1-8-7) are ablation artifacts only.
    // They are not in FORMAT-SPEC-001 and must not be confused with the
    // canonical format. The ablation is additive and read-only w.r.t. L6.
    @compileAssert(ALT_A_EXP_BITS != GF16_CANONICAL_EXP_BITS);
    @compileAssert(ALT_A_MANT_BITS != GF16_CANONICAL_MANT_BITS);
}

invariant "gf16_lp_canonical_minimises_phi_dist" {
    // [Registered structural fact from FORMAT-SPEC-001]
    // phi_dist(1-6-9) < phi_dist(1-7-8); the canonical split is closer to phi^-1.
    // This is a STRUCTURAL fact, not a training quality claim.
    @compileAssert(GF16_PHI_DIST_REGISTERED < ALT_A_PHI_DIST);
}

// ============================================================================
// C3 Bench (L4 TESTABILITY)
// ============================================================================

bench "bench_gf16_lp_encode_fp8_projection" {
    // Latency of: encode f32 to GF16, truncate mantissa to 4 bits (FP8 floor),
    //   decode back. This is the inner loop cost of FP8-projected GF16 operations.
    // Target: < 2x overhead vs plain gf16_encode_f32 + gf16_decode_to_f32.
    @setEvalBranchQuota(10000);
    const v : f32 = 1.5;
    const fp8_shift : u4 = GF16_CANONICAL_MANT_BITS - FP8_EFFECTIVE_MANT_BITS;
    var result : f32 = 0.0;
    for (0..1000) |_| {
        const enc = gf16_encode_f32(v);
        const mant = enc & MANT_MASK;
        const mant_trunc = (mant >> fp8_shift) << fp8_shift;
        const enc_trunc = (enc & ~MANT_MASK) | mant_trunc;
        result = gf16_decode_to_f32(enc_trunc);
    }
    _ = result;
}

bench "bench_gf16_lp_roundtrip_baseline" {
    // Plain GF16 encode + decode (no FP8 truncation) for latency comparison.
    // The ratio bench_gf16_lp_encode_fp8_projection / bench_gf16_lp_roundtrip_baseline
    //   must be < 2.0 (the MAX_LATENCY_RATIO from C2 convention).
    @setEvalBranchQuota(10000);
    const v : f32 = 1.5;
    var result : f32 = 0.0;
    for (0..1000) |_| {
        const enc = gf16_encode_f32(v);
        result = gf16_decode_to_f32(enc);
    }
    _ = result;
}

Open the lesson's spec in the player ↗

All lessons

• Module 1 · The rule and its numbers

  One rule splits every width, the ratio it aims at, and the Lucas numbers behind the 3.

  1. 01 A float cut by phi
  2. 02 Phi as a ratio
  3. 03 Lucas numbers stay whole

• Module 2 · Why phi, why three

  Why the split is phi, why base three, and how a spec checks GF16 keeps phi.

  1. 04 Why the split is phi
  2. 05 Why base three
  3. 06 Phi in sixteen bits

• Module 3 · The small rungs: GF4 to GF8

  GF4, GF6 and GF8, the fewest bits, where rounding to whole bits costs the most.

  1. 07 GF4: four bits
  2. 08 GF6: six bits
  3. 09 GF8: one byte

• Module 4 · Ten to fourteen bits

  GF10, GF12 and GF14, and how the distance from 1 / phi moves as the word grows.

  1. 10 GF10: ten bits
  2. 11 GF12: twelve bits
  3. 12 GF14: fourteen bits

• Module 5 · GF16 at work

  The primary 16-bit format, a two-term dot product in GF-T16, then GF20 and GF24.

  1. 13 GF16: the primary format
  2. 14 Two products, one sum
  3. 15 GF20 and GF24

• Module 6 · GF32 to GF64

  GF32 beside IEEE single, GF48 with no IEEE twin, GF64 beside IEEE double.

  1. 16 GF32: a single
  2. 17 GF48: forty-eight bits
  3. 18 GF64: a double

• Module 7 · GF96 to GF256

  GF96, GF128 and GF256, where the specs hold the layout with invariants.

  1. 19 GF96: ninety-six bits
  2. 20 GF128: a quad
  3. 21 GF256: two hundred fifty-six bits

• Module 8 · The widest rungs, then trits

  GF512 and GF1024, the two widest rungs, then GF-T8, where the exponent moves to trits.

  1. 22 GF512: five hundred twelve bits
  2. 23 GF1024: one kilobit
  3. 24 GF-T8: an exponent in trits

• Module 9 · More trits, then the decode

  GF-T16 and GF-T32, then why fixed fields decode in parallel and a posit does not.

  1. 25 GF-T16: six trits
  2. 26 GF-T32: twelve trits
  3. 27 Decode in one step