//! Q16_16 fixed-point arithmetic — Rust Port
//!
//! Mirrors `Core/SilverSight/FixedPoint.lean` (Lean 4).
//!
//! Rounding mode: Banker's rounding (round-half-to-even), matching Python `round()`.
//! This differs from Lean `Int.div` (floor division). The discrepancy is documented
//! and cross-validated against Python, R, and Julia outputs.
//!
//! Design rules:
//!   - Core type is `i32` (signed 32-bit, matches `Q16_16Value`)
//!   - Saturation via `Self::of_raw`
//!   - Float is external boundary only (`of_float`, `to_float`)

use std::cmp::{max, min};

pub const Q16_MIN_RAW: i32 = -2147483648;
pub const Q16_MAX_RAW: i32 = 2147483647;
pub const Q16_SCALE: i32 = 65536;

/// Q16.16 fixed-point raw integer, saturated into `[Q16_MIN_RAW, Q16_MAX_RAW]`.
pub type Q16 = i32;

/// Saturating clamp matching Lean `q16Clamp`.
#[inline]
pub fn clamp_raw(raw: i64) -> i32 {
    if raw > Q16_MAX_RAW as i64 {
        Q16_MAX_RAW
    } else if raw < Q16_MIN_RAW as i64 {
        Q16_MIN_RAW
    } else {
        raw as i32
    }
}

/// Saturating constructor. Matches Lean `Q16_16.ofRawInt`.
#[inline]
pub fn of_raw(raw: i32) -> Q16 {
    raw // already in-range if constructed via clamp_raw
}

/// Float → Q16_16 (external boundary only). Matches Python `f2q` (banker's rounding).
pub fn of_float(f: f64) -> Q16 {
    if f.is_nan() || f >= 32768.0 {
        return Q16_MAX_RAW;
    }
    if f <= -32768.0 {
        return Q16_MIN_RAW;
    }
    clamp_raw((f * 65536.0).round_ties_even() as i64)
}

/// Q16_16 → Float64.
pub fn to_float(q: Q16) -> f64 {
    q as f64 / 65536.0
}

/// Constants (match Lean `#eval` witnesses)
pub const ZERO: Q16 = 0;
pub const ONE: Q16 = 65536;
pub const NEG_ONE: Q16 = -65536;
pub const EPSILON: Q16 = 1;
pub const TWO: Q16 = 131072;
pub const HALF: Q16 = 32768;
pub const MAX_VAL: Q16 = Q16_MAX_RAW;
pub const MIN_VAL: Q16 = Q16_MIN_RAW;

// NUVMAP constants
pub const PCT1: Q16 = 655;      // 0.01 floor
pub const PCT70: Q16 = 45875;   // 0.7 floor
pub const PCT30: Q16 = 19661;   // 0.3 ceil (banker's)
pub const ONE50: Q16 = 98304;   // 1.5 exact

// Mathematical constants (from Q16_16Numerics.lean)
pub const PI: Q16 = 205887;
pub const E: Q16 = 178145;
pub const LN2: Q16 = 45426;

// ── Arithmetic (saturating, banker's rounding) ─────────────────────────

/// Saturating addition. Matches Lean `Q16_16.add`.
#[inline]
pub fn add(a: Q16, b: Q16) -> Q16 {
    a.saturating_add(b)
}

/// Saturating subtraction.
#[inline]
pub fn sub(a: Q16, b: Q16) -> Q16 {
    a.saturating_sub(b)
}

/// Q16.16 multiplication: `round((a * b) / 65536)`. Uses i64 intermediate.
pub fn mul(a: Q16, b: Q16) -> Q16 {
    let prod = (a as i64) * (b as i64);
    // Banker's rounding: round half to even
    let scaled = ((prod as f64) / (Q16_SCALE as f64)).round_ties_even() as i64;
    clamp_raw(scaled)
}

/// Q16.16 division: `round((a * 65536) / b)`. Returns MAX_VAL if b == 0.
pub fn div_op(a: Q16, b: Q16) -> Q16 {
    if b == 0 {
        return MAX_VAL;
    }
    let prod = (a as i64) * (Q16_SCALE as i64);
    let scaled = ((prod as f64) / (b as f64)).round_ties_even() as i64;
    clamp_raw(scaled)
}

/// Saturating negation.
#[inline]
pub fn neg(a: Q16) -> Q16 {
    a.saturating_neg()
}

/// Absolute value.
#[inline]
pub fn abs_(a: Q16) -> Q16 {
    if a < 0 { neg(a) } else { a }
}

// ── Comparisons ────────────────────────────────────────────────────────

#[inline]
pub fn lt(a: Q16, b: Q16) -> bool { a < b }
#[inline]
pub fn le(a: Q16, b: Q16) -> bool { a <= b }
#[inline]
pub fn gt(a: Q16, b: Q16) -> bool { a > b }
#[inline]
pub fn ge(a: Q16, b: Q16) -> bool { a >= b }

// ── Numeric functions ───────────────────────────────────────────────────

/// Integer square root via Newton's method (floor). Matches Lean `intSqrt`.
pub fn int_sqrt(n: i64) -> i64 {
    if n <= 0 { return 0; }
    let mut x = n / 2 + 1;
    for _ in 0..64 {
        let x_new = (x + n / x) / 2;
        if x_new >= x { return x; }
        x = x_new;
    }
    x
}

/// Q16.16 square root.
pub fn sqrt(q: Q16) -> Q16 {
    if q <= 0 { return ZERO; }
    clamp_raw(int_sqrt((q as i64) * (Q16_SCALE as i64)) as i64)
}

/// Natural logarithm approximation around 1.0. Matches Lean `Q16_16.ln`.
pub fn ln_nat(q: Q16) -> Q16 {
    let x = q as i64;
    if x <= 0 { return ZERO; }
    let y = x - Q16_SCALE as i64;
    let y2 = (y * y) / Q16_SCALE as i64;
    let y3 = (y * y2) / Q16_SCALE as i64;
    clamp_raw((y - y2 / 2 + y3 / 3) as i64)
}

/// Q16.16 exponential via Taylor with range reduction. Matches Lean `Q16_16.exp`.
pub fn exp(x: Q16) -> Q16 {
    let rx = x as i64;
    if rx <= -(Q16_SCALE as i64) { return ZERO; }
    if rx >= 4 * (Q16_SCALE as i64) { return MAX_VAL; }

    let ln2_raw: i64 = 45426;
    let k = rx / ln2_raw;
    let r = rx - k * ln2_raw;
    let scale = Q16_SCALE as i64;

    let r2 = (r * r) / scale;
    let r3 = (r * r2) / scale;
    let r4 = (r2 * r2) / scale;
    let r5 = (r2 * r3) / scale;
    let r6 = (r3 * r3) / scale;

    let taylor = scale + r + r2 / 2 + r3 / 6 + r4 / 24 + r5 / 120 + r6 / 720;

    let shifted = if k >= 15 { Q16_MAX_RAW as i64 }
    else if k <= -15 { 0 }
    else if k >= 0 { taylor << (k as u32) }
    else { taylor >> ((-k) as u32) };

    clamp_raw(shifted)
}

// ── Tests ───────────────────────────────────────────────────────────────

#[cfg(test)]
mod tests {
    use crate::q16::*;

    #[test]
    fn test_constants() {
        assert_eq!(ZERO, 0);
        assert_eq!(ONE, 65536);
        assert_eq!(EPSILON, 1);
        assert_eq!(HALF, 32768);
        assert_eq!(to_float(PCT1), 0.0099945068359375);
        assert_eq!(to_float(PCT70), 0.6999969482421875);
        assert_eq!(to_float(PCT30), 0.3000030517578125);
        assert_eq!(to_float(ONE50), 1.5);
    }

    #[test]
    fn test_arithmetic() {
        // 1.0 * 1.0 = 1.0
        let r = mul(ONE, ONE);
        assert_eq!(to_float(r), 1.0);

        // 1.0 / 2.0 = 0.5
        let r = div_op(ONE, TWO);
        assert_eq!(to_float(r), 0.5);

        // 0.5 + 0.5 = 1.0
        let r = add(HALF, HALF);
        assert_eq!(r, ONE);

        // 1.0 - 0.5 = 0.5
        let r = sub(ONE, HALF);
        assert_eq!(r, HALF);
    }

    #[test]
    fn test_of_float() {
        // of_float and to_float round-trip for simple cases
        let q = of_float(1.0);
        assert_eq!(q, ONE);
        assert!((to_float(q) - 1.0).abs() < 1e-6);

        let q = of_float(0.5);
        assert_eq!(q, HALF);

        // f2q(0.8) = round(52428.8) = 52429
        let q = of_float(0.8);
        assert_eq!(q, 52429);

        // f2q(0.75) = round(49152.0) = 49152
        let q = of_float(0.75);
        assert_eq!(q, 49152);
    }

    #[test]
    fn test_nuvmap_data_baseline() {
        // Cross-validated against Python/R/Julia for cell 0:
        // amvr=0.8, avmr=0.75, achiral_stable
        let amvr_q = of_float(0.8);
        let avmr_q = of_float(0.75);
        assert_eq!(amvr_q, 52429);
        assert_eq!(avmr_q, 49152);

        // raw_eigenmass = q_div(101581, 131072) = 50790
        let raw = div_op(add(amvr_q, avmr_q), of_raw(131072));
        assert_eq!(raw, 50790);

        // max_eigenmass (from cell 4)
        let amvr4 = of_float(0.9);
        let avmr4 = of_float(0.85);
        let max_raw = div_op(add(amvr4, avmr4), of_raw(131072));
        assert_eq!(max_raw, 57344);

        // E_norm
        let e_norm = div_op(raw, max_raw);
        assert_eq!(e_norm, 58046);
        assert!((to_float(e_norm) - 0.885711669921875).abs() < 1e-10);

        // R_i
        let one_minus = sub(ONE, e_norm);
        let r_i = if one_minus > PCT1 { one_minus } else { PCT1 };
        assert_eq!(r_i, 7490);

        // E_i
        let e_i = div_op(mul(mul(mul(ONE, e_norm), ONE), ONE), add(r_i, EPSILON));
        assert_eq!(e_i, 507823);
        assert!((to_float(e_i) - 7.7487640380859375).abs() < 1e-10);
    }
}
