g_math 0.6.6

//! The f64/f32 boundary of `FixedPoint`.
//!
//! `to_f64` / `to_f32` are exact when the raw value fits the format's
//! significand and rounded to nearest with ties to even otherwise (f32:
//! subnormal or zero below its normal range, infinite above it).
//! `from_f64` / `from_f32` truncate toward zero; `try_from_f64` /
//! `try_from_f32` return `InvalidInput` for NaN and `TierOverflow` for infinity
//! or a value outside the profile's range, and the infallible forms panic there
//! instead of wrapping.
//!
//! References: `tests/data/float_boundary_refs.rs`, generated by
//! `scripts/generate_float_boundary_refs.py` from exact rationals (round half
//! to even), f64 cross-checked against CPython's correctly rounded division.
//! Every comparison here is on bit patterns and raw integers; no float
//! arithmetic is performed.

use g_math::fixed_point::{FixedPoint, OverflowDetected};

#[allow(dead_code)]
mod refs {
    include!("data/float_boundary_refs.rs");
}
use refs::{FromFloat, ToFloat};

#[cfg(table_format = "q128_128")]
use g_math::fixed_point::domains::binary_fixed::i256::I256;
#[cfg(table_format = "q256_256")]
use g_math::fixed_point::domains::binary_fixed::i512::I512;

// ============================================================================
// Raw values from the references, and the profile's reference tables
// ============================================================================

/// Little-endian 64-bit words of a hex magnitude.
#[cfg(any(table_format = "q128_128", table_format = "q256_256"))]
fn hex_words(hex: &str, count: usize) -> Vec<u64> {
    let mut words = vec![0u64; count];
    let digits = hex.as_bytes();
    for (i, chunk) in digits.rchunks(16).enumerate() {
        words[i] = u64::from_str_radix(std::str::from_utf8(chunk).unwrap(), 16).unwrap();
    }
    words
}

#[cfg(any(table_format = "q16_16", table_format = "q32_32", table_format = "q64_64"))]
fn raw_from(negative: bool, hex: &str) -> FixedPoint {
    let magnitude = u128::from_str_radix(hex, 16).unwrap();
    let value = if negative { (magnitude as i128).wrapping_neg() } else { magnitude as i128 };
    FixedPoint::from_raw(value as _)
}

#[cfg(table_format = "q128_128")]
fn raw_from(negative: bool, hex: &str) -> FixedPoint {
    let words = hex_words(hex, 4);
    let magnitude = I256 { words: [words[0], words[1], words[2], words[3]] };
    FixedPoint::from_raw(if negative { -magnitude } else { magnitude })
}

#[cfg(table_format = "q256_256")]
fn raw_from(negative: bool, hex: &str) -> FixedPoint {
    let words = hex_words(hex, 8);
    let mut array = [0u64; 8];
    array.copy_from_slice(&words);
    let magnitude = I512 { words: array };
    FixedPoint::from_raw(if negative { -magnitude } else { magnitude })
}

/// A small raw value (|raw| < 2^63) on any profile.
fn small_raw(value: i64) -> FixedPoint {
    #[cfg(any(table_format = "q16_16", table_format = "q32_32", table_format = "q64_64"))]
    {
        FixedPoint::from_raw(value as _)
    }
    #[cfg(any(table_format = "q128_128", table_format = "q256_256"))]
    {
        raw_from(value < 0, &format!("{:x}", value.unsigned_abs()))
    }
}

/// Fraction bits of this build: one is `2^F` raw, so count the exact halvings
/// that take it down to raw 1.
fn frac_bits() -> u32 {
    // the build's split (halving one with from_int(2) fails where 2 is out of
    // range, at GMATH_FRAC_BITS = 30)
    g_math::fixed_point::frac_config::FRAC_BITS
}

/// The profile's reference tables, or `None` for a realtime split the
/// generator does not cover.
fn tables() -> Option<(&'static [ToFloat], &'static [FromFloat], &'static [FromFloat])> {
    #[cfg(table_format = "q16_16")]
    {
        match frac_bits() {
            16 => Some((refs::REALTIME_F16_TO_FLOAT, refs::REALTIME_F16_FROM_F64, refs::REALTIME_F16_FROM_F32)),
            10 => Some((refs::REALTIME_F10_TO_FLOAT, refs::REALTIME_F10_FROM_F64, refs::REALTIME_F10_FROM_F32)),
            other => {
                println!("no references for GMATH_FRAC_BITS={other}: add it to the generator");
                None
            }
        }
    }
    #[cfg(table_format = "q32_32")]
    {
        Some((refs::COMPACT_TO_FLOAT, refs::COMPACT_FROM_F64, refs::COMPACT_FROM_F32))
    }
    #[cfg(table_format = "q64_64")]
    {
        Some((refs::EMBEDDED_TO_FLOAT, refs::EMBEDDED_FROM_F64, refs::EMBEDDED_FROM_F32))
    }
    #[cfg(table_format = "q128_128")]
    {
        Some((refs::BALANCED_TO_FLOAT, refs::BALANCED_FROM_F64, refs::BALANCED_FROM_F32))
    }
    #[cfg(table_format = "q256_256")]
    {
        Some((refs::SCIENTIFIC_TO_FLOAT, refs::SCIENTIFIC_FROM_F64, refs::SCIENTIFIC_FROM_F32))
    }
}

// ============================================================================
// Reference gates
// ============================================================================

#[test]
fn to_float_matches_the_exact_rational_references() {
    let Some((to_float, _, _)) = tables() else { return };
    for r in to_float {
        let x = raw_from(r.negative, r.magnitude);
        let sign = if r.negative { "-" } else { "" };
        assert_eq!(x.to_f64().to_bits(), r.f64_bits, "to_f64 of raw {sign}0x{}", r.magnitude);
        assert_eq!(x.to_f32().to_bits(), r.f32_bits, "to_f32 of raw {sign}0x{}", r.magnitude);
    }
    println!("{} raws: to_f64 and to_f32 bit-identical to the references", to_float.len());
}

fn check_from(kind: &str, vectors: &[FromFloat], convert: impl Fn(u64) -> Result<FixedPoint, OverflowDetected>) {
    for v in vectors {
        match (convert(v.bits), v.raw) {
            (Ok(x), Some((negative, hex))) => {
                assert!(x == raw_from(negative, hex), "{kind} of bits 0x{:x}: raw differs from the reference", v.bits)
            }
            (Err(OverflowDetected::TierOverflow), None) => {}
            (got, expected) => panic!(
                "{kind} of bits 0x{:x}: got {:?}, reference {}",
                v.bits,
                got.map(|x| x.raw()),
                if expected.is_some() { "a raw value" } else { "TierOverflow" }
            ),
        }
    }
    println!("{kind}: {} inputs match the truncation references", vectors.len());
}

#[test]
fn from_float_truncates_toward_zero_and_rejects_out_of_range() {
    let Some((_, from_f64, from_f32)) = tables() else { return };
    check_from("try_from_f64", from_f64, |bits| FixedPoint::try_from_f64(f64::from_bits(bits)));
    check_from("try_from_f32", from_f32, |bits| FixedPoint::try_from_f32(f32::from_bits(bits as u32)));
    // the infallible forms agree wherever the conversion succeeds
    for v in from_f64.iter().filter(|v| v.raw.is_some()) {
        let x = f64::from_bits(v.bits);
        assert!(FixedPoint::from_f64(x) == FixedPoint::try_from_f64(x).unwrap());
    }
}

// ============================================================================
// Exactness and round trips (no references needed)
// ============================================================================

/// `bits` (an f64) equals `raw x 2^-frac` exactly, checked on integers.
fn f64_is_exactly(bits: u64, raw: i64, frac: u32) -> bool {
    let negative = bits >> 63 == 1;
    let biased = ((bits >> 52) & 0x7FF) as i32;
    if raw == 0 {
        return bits == 0;
    }
    if negative != (raw < 0) || biased == 0 || biased == 0x7FF {
        return false;
    }
    let significand = (bits & ((1u64 << 52) - 1)) | (1u64 << 52);
    // value = significand x 2^(biased - 1075); compare with |raw| x 2^-frac
    let shift = biased - 1075 + frac as i32;
    let magnitude = raw.unsigned_abs() as u128;
    if shift >= 0 {
        (significand as u128).checked_shl(shift as u32) == Some(magnitude)
    } else {
        magnitude.checked_shl((-shift) as u32) == Some(significand as u128)
    }
}

/// Deterministic raws with |raw| < 2^53 (and within the profile's range).
fn sample_raws() -> Vec<i64> {
    #[cfg(table_format = "q16_16")]
    let bound: i64 = 1 << 31;
    #[cfg(not(table_format = "q16_16"))]
    let bound: i64 = 1 << 53;
    let mut raws: Vec<i64> = (-65_536..=65_536).collect();
    let mut state = 0x2026_0919_u64;
    for _ in 0..100_000 {
        state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
        let bits = 1 + (state >> 58) as u32 % 53; // magnitude bit length 1..=53
        state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
        let magnitude = ((state >> 11) as i64) & ((1i64 << bits) - 1);
        let value = if state & 1 == 1 { -magnitude } else { magnitude };
        if value > -bound && value < bound {
            raws.push(value);
        }
    }
    raws
}

#[test]
fn exact_conversions_are_exact_and_round_trip() {
    let frac = frac_bits();
    let mut f32_checked = 0;
    for raw in sample_raws() {
        let x = small_raw(raw);
        let y = x.to_f64();
        assert!(f64_is_exactly(y.to_bits(), raw, frac), "to_f64 of raw {raw} is not exact");
        assert!(FixedPoint::from_f64(y) == x, "from_f64(to_f64(x)) != x for raw {raw}");
        // exact in f32 when it has at most 24 bits and is a multiple of f32's
        // smallest subnormal 2^-149 (not on scientific, whose step is 2^-256)
        if raw.unsigned_abs() < 1 << 24 && frac <= 149 {
            assert!(FixedPoint::from_f32(x.to_f32()) == x, "from_f32(to_f32(x)) != x for raw {raw}");
            f32_checked += 1;
        }
    }
    println!("round trips exact: {} f64, {f32_checked} f32", sample_raws().len());
}

#[test]
fn round_trips_do_not_erode() {
    // 0.6.2 lost one ulp toward zero per round trip (99, 98, 97, ... at Q22.10)
    for raw in [99i64, -99, 1, -1, 12_345, -12_345] {
        let x = small_raw(raw);
        let mut y = FixedPoint::from_f64(x.to_f64());
        for _ in 0..5 {
            y = FixedPoint::from_f64(y.to_f64());
        }
        assert!(y == x, "raw {raw} moved after repeated round trips");
    }
}

#[test]
fn to_f64_is_sign_symmetric_and_monotone() {
    let mut raws: Vec<i64> = sample_raws().into_iter().filter(|r| *r > 0).collect();
    raws.sort_unstable();
    raws.dedup();
    let mut previous = 0u64;
    for raw in raws {
        let positive = small_raw(raw).to_f64().to_bits();
        let negative = small_raw(-raw).to_f64().to_bits();
        assert_eq!(negative, positive | (1 << 63), "to_f64(-x) != -to_f64(x) for raw {raw}");
        // positive f64 bit patterns order like their values
        assert!(positive > previous, "to_f64 not strictly increasing at raw {raw}");
        previous = positive;
    }
}

// ============================================================================
// Errors and panics
// ============================================================================

#[test]
fn non_finite_inputs_are_errors() {
    assert_eq!(FixedPoint::try_from_f64(f64::NAN).unwrap_err(), OverflowDetected::InvalidInput);
    assert_eq!(FixedPoint::try_from_f32(f32::NAN).unwrap_err(), OverflowDetected::InvalidInput);
    for v in [f64::INFINITY, f64::NEG_INFINITY] {
        assert_eq!(FixedPoint::try_from_f64(v).unwrap_err(), OverflowDetected::TierOverflow);
    }
    for v in [f32::INFINITY, f32::NEG_INFINITY] {
        assert_eq!(FixedPoint::try_from_f32(v).unwrap_err(), OverflowDetected::TierOverflow);
    }
    assert!(FixedPoint::try_from_f64(-0.0).unwrap() == FixedPoint::ZERO);
}

/// `2^e` as f64 bits (for -1022 <= e <= 1023).
fn power_of_two(e: i32) -> f64 {
    f64::from_bits(((e + 1023) as u64) << 52)
}

/// The profile's range is `[-2^top, 2^top)`.
fn range_top() -> i32 {
    #[cfg(table_format = "q16_16")]
    let width = 32;
    #[cfg(table_format = "q32_32")]
    let width = 64;
    #[cfg(table_format = "q64_64")]
    let width = 128;
    #[cfg(table_format = "q128_128")]
    let width = 256;
    #[cfg(table_format = "q256_256")]
    let width = 512;
    width - 1 - frac_bits() as i32
}

#[test]
fn the_range_ends_exactly_at_its_bounds() {
    let top = range_top();
    // -2^top is the storage minimum; 2^top is one step beyond the maximum
    let minimum = FixedPoint::try_from_f64(-power_of_two(top)).expect("-2^top is representable");
    assert!(minimum < FixedPoint::ZERO);
    assert_eq!(FixedPoint::try_from_f64(power_of_two(top)).unwrap_err(), OverflowDetected::TierOverflow);
    assert_eq!(FixedPoint::try_from_f64(-power_of_two(top + 1)).unwrap_err(), OverflowDetected::TierOverflow);
}

#[test]
#[should_panic(expected = "infinity or value outside the profile's range")]
fn from_f64_panics_beyond_the_range_instead_of_wrapping() {
    FixedPoint::from_f64(power_of_two(range_top()));
}

#[test]
#[should_panic(expected = "NaN")]
fn from_f64_panics_on_nan() {
    FixedPoint::from_f64(f64::NAN);
}

#[test]
#[should_panic(expected = "infinity or value outside the profile's range")]
fn from_f32_panics_on_infinity() {
    FixedPoint::from_f32(f32::INFINITY);
}