solmath 0.2.0

#![allow(dead_code)]

// ============================================================
// Fixed-point constants
// ============================================================

/// Standard fixed-point scale factor: 1e12.
///
/// All unsigned fixed-point values are integers representing `value * SCALE`.
pub const SCALE: u128 = 1_000_000_000_000u128; // 1e12

/// Signed standard scale factor: 1e12.
///
/// Same magnitude as [`SCALE`] but typed as `i128` for signed arithmetic.
pub const SCALE_I: i128 = 1_000_000_000_000i128; // 1e12
pub const SCALE_6: i64 = 1_000_000i64; // 1e6
pub const SCALE_HP_U: u128 = 1_000_000_000_000_000u128; // 1e15
pub const SCALE_HP: i128 = 1_000_000_000_000_000i128; // 1e15
pub const HP_TO_STD: i128 = 1_000i128; // 1e15 / 1e12
pub const LN2: u128 = 693_147_180_560u128; // ln(2) * 1e12
pub const LN2_I: i128 = 693_147_180_560i128; // ln(2) * 1e12
pub const LN2_HP: i128 = 693_147_180_559_945i128; // ln(2) * 1e15
pub const HALF_LN2_HP: i128 = 346_573_590_279_973i128; // ln(2)/2 * 1e15
pub const INV_SQRT2: u128 = 707_106_781_187u128; // 1/sqrt(2) * 1e12
pub const SQRT2_SCALE: u128 = 1_414_213_562_373u128; // sqrt(2) * 1e12
pub const INV_SQRT2_HP: u128 = 707_106_781_186_548u128; // 1/sqrt(2) * 1e15
pub const SQRT2_SCALE_HP: u128 = 1_414_213_562_373_095u128; // sqrt(2) * 1e15
pub const INV_SQRT_2PI_HP: i128 = 398_942_280_401_433i128; // 1/sqrt(2pi) * 1e15
pub const LN_2PI: i128 = 1_837_877_066_409i128; // ln(2π) * 1e12
pub const LN_HP_TERMS: usize = 12;
pub const EXP_HP_TERMS: i128 = 15;
pub const LN2_6: i64 = 693_147; // ln(2) * 1e6
pub const PI6: i64 = 3_141_593;
pub const PI2_6: i64 = 6_283_185;
pub const PIH_6: i64 = 1_570_796;
pub const PIQ_6: i64 = 785_398;
pub const SC0_6: i64 = 1_000_000;
pub const SC1_6: i64 = -166_667;
pub const SC2_6: i64 = 8_333;
pub const SC3_6: i64 = -198;
pub const SC4_6: i64 = 3;
pub const CC0_6: i64 = 1_000_000;
pub const CC1_6: i64 = -500_000;
pub const CC2_6: i64 = 41_667;
pub const CC3_6: i64 = -1_389;
pub const CC4_6: i64 = 25;
pub const FP_DIV_THIN_MAX: u128 = u128::MAX / SCALE;

// Remez minimax polynomial for ln via arctanh substitution (wide domain).
// P(u) ~ arctanh(sqrt(u))/sqrt(u) on u in [0, 1/9], degree 7.
// Usage: ln(m/SCALE) = 2*t*P(t^2) where t = (m-SCALE)/(m+SCALE).
// No sqrt(2) secondary reduction needed — covers full [SCALE, 2*SCALE).
// Approximation error: < 0.05 ULP. Generated by scripts/fit_ln_remez_wide.py.
pub const LN_REMEZ_W0: i128 = 1_000_000_000_000;
pub const LN_REMEZ_W1: i128 = 333_333_333_406;
pub const LN_REMEZ_W2: i128 = 199_999_986_321;
pub const LN_REMEZ_W3: i128 = 142_858_120_138;
pub const LN_REMEZ_W4: i128 = 111_077_000_901;
pub const LN_REMEZ_W5: i128 = 91_550_218_194;
pub const LN_REMEZ_W6: i128 = 70_353_796_623;
pub const LN_REMEZ_W7: i128 = 100_340_778_474;

// HP Remez minimax polynomial for ln at SCALE_HP = 1e15 (wide domain).
// Degree 9, approximation error < 0.033 ULP. Generated by scripts/fit_ln_remez_wide_hp.py.
pub const LN_REMEZ_HP0: i128 = 1_000_000_000_000_000;
pub const LN_REMEZ_HP1: i128 = 333_333_333_333_413;
pub const LN_REMEZ_HP2: i128 = 199_999_999_976_438;
pub const LN_REMEZ_HP3: i128 = 142_857_145_558_286;
pub const LN_REMEZ_HP4: i128 = 111_110_954_306_731;
pub const LN_REMEZ_HP5: i128 = 90_914_299_002_370;
pub const LN_REMEZ_HP6: i128 = 76_818_381_244_458;
pub const LN_REMEZ_HP7: i128 = 67_956_527_427_899;
pub const LN_REMEZ_HP8: i128 = 49_403_945_243_824;
pub const LN_REMEZ_HP9: i128 = 88_903_819_130_567;

/// Sub-ULP residual for split LN2 at standard scale.
/// true_ln2 × SCALE ≈ LN2_I + LN2_LO / SCALE to sub-ULP precision.
/// Replaces the crude -(k*55+500)/1000 correction in ln and exp.
pub const LN2_LO: i128 = -54_690_582_768;

/// Sub-ULP residual for split LN2 at HP scale.
pub const LN2_HP_LO: i128 = 309_417_232_121_458;

/// Remez degree-7 ln polynomial as array for compensated evaluation.
/// Same coefficients as LN_REMEZ_W0..W7, ascending order.
pub const LN_REMEZ_COEFFS: [i128; 8] = [
    LN_REMEZ_W0,
    LN_REMEZ_W1,
    LN_REMEZ_W2,
    LN_REMEZ_W3,
    LN_REMEZ_W4,
    LN_REMEZ_W5,
    LN_REMEZ_W6,
    LN_REMEZ_W7,
];

/// HP Remez degree-9 ln polynomial as array for compensated evaluation.
/// Same coefficients as LN_REMEZ_HP0..HP9, ascending order.
pub const LN_REMEZ_HP_COEFFS: [i128; 10] = [
    LN_REMEZ_HP0,
    LN_REMEZ_HP1,
    LN_REMEZ_HP2,
    LN_REMEZ_HP3,
    LN_REMEZ_HP4,
    LN_REMEZ_HP5,
    LN_REMEZ_HP6,
    LN_REMEZ_HP7,
    LN_REMEZ_HP8,
    LN_REMEZ_HP9,
];

// Remez rational coefficients for exp.
// c(r) = r - P(r^2)*r^2, exp(r) = 1 + r + r*c/(2-c).
// From FreeBSD msun, degree 5 in r^2. Generated by scripts/fit_exp_remez.py.
pub const EXP_REMEZ_P1: i128 = 166_666_666_667;
pub const EXP_REMEZ_P2: i128 = -2_777_777_778;
pub const EXP_REMEZ_P3: i128 = 66_137_563;
pub const EXP_REMEZ_P4: i128 = -1_653_390;
pub const EXP_REMEZ_P5: i128 = 41_381;

// HP Remez rational coefficients for exp at SCALE_HP = 1e15.
pub const EXP_REMEZ_HP_P1: i128 = 166_666_666_666_666;
pub const EXP_REMEZ_HP_P2: i128 = -2_777_777_777_702;
pub const EXP_REMEZ_HP_P3: i128 = 66_137_563_214;
pub const EXP_REMEZ_HP_P4: i128 = -1_653_390_221;
pub const EXP_REMEZ_HP_P5: i128 = 41_381_368;

// ============================================================
// MulDivRounding enum
// ============================================================

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum MulDivRounding {
    ToZero,
    Floor,
    Ceil,
}

// ============================================================
// U256 struct and impl
// ============================================================

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct U256 {
    pub limbs: [u64; 4],
}

/// Assemble a 256-bit product from the four exact 64×64 partial products.
///
/// Keeping carry propagation separate from partial-product generation gives
/// the formal harness a solver-sized boundary: the harness can quantify over
/// every valid partial product and prove each output limb independently.
#[inline]
fn assemble_u128_product(p00: u128, p01: u128, p10: u128, p11: u128) -> U256 {
    const MASK: u128 = u64::MAX as u128;

    // Assemble one base-2^64 column at a time. Every valid column is below
    // 3*2^64, so it fits comfortably in u128 even though the unsplit middle
    // partial products would not fit when added directly.
    let column1 = (p00 >> 64) + (p01 & MASK) + (p10 & MASK);
    let column2 = (p01 >> 64) + (p10 >> 64) + (p11 & MASK) + (column1 >> 64);
    let column3 = (p11 >> 64) + (column2 >> 64);

    U256 {
        limbs: [p00 as u64, column1 as u64, column2 as u64, column3 as u64],
    }
}

/// One base-2^64 long-division step.
///
/// The precondition `remainder < divisor` makes the quotient digit fit in one
/// limb. `U256::div_rem_u64` applies this exact transition four times.
#[inline]
fn div_rem_u64_step(remainder: u64, limb: u64, divisor: u64) -> (u64, u64) {
    debug_assert!(divisor != 0);
    debug_assert!(remainder < divisor);
    let current = ((remainder as u128) << 64) | limb as u128;
    (
        (current / divisor as u128) as u64,
        (current % divisor as u128) as u64,
    )
}

/// Knuth Algorithm D quotient-digit refinement for a normalized two-limb
/// divisor. Normalization guarantees `v1 >= 2^63`, so at most two decrements
/// are required before the trial digit is valid or `rhat` crosses the base.
#[inline]
fn refine_knuth_quotient_digit(
    mut qhat: u128,
    mut rhat: u128,
    u0: u64,
    v0: u64,
    v1: u64,
) -> (u128, u128) {
    const BASE: u128 = 1u128 << 64;
    debug_assert!(qhat < BASE);
    debug_assert!(u128::from(v1) >= BASE / 2);

    if rhat < BASE {
        while qhat * u128::from(v0) > (rhat << 64) + u128::from(u0) {
            debug_assert!(qhat > 0);
            (qhat, rhat) = advance_knuth_refinement(qhat, rhat, v1);
            if rhat >= BASE {
                break;
            }
        }
    }
    (qhat, rhat)
}

/// Apply the arithmetic state transition for one Knuth-D3 correction.
#[inline]
fn advance_knuth_refinement(qhat: u128, rhat: u128, v1: u64) -> (u128, u128) {
    const BASE: u128 = 1u128 << 64;
    debug_assert!(qhat > 0);
    debug_assert!(rhat < BASE);
    debug_assert!(u128::from(v1) >= BASE / 2);
    (qhat - 1, rhat + u128::from(v1))
}

impl U256 {
    #[inline]
    pub const fn zero() -> Self {
        Self { limbs: [0; 4] }
    }

    /// Create U256 from a u128 value.
    #[inline]
    pub fn from_u128(x: u128) -> Self {
        Self {
            limbs: [x as u64, (x >> 64) as u64, 0, 0],
        }
    }

    /// Extract the low 128 bits.
    #[inline]
    pub fn low_u128(&self) -> u128 {
        (self.limbs[0] as u128) | ((self.limbs[1] as u128) << 64)
    }

    /// Check if the high 128 bits are nonzero.
    #[inline]
    pub fn high_u128_nonzero(&self) -> bool {
        self.limbs[2] != 0 || self.limbs[3] != 0
    }

    /// Get bit at index.
    #[inline]
    pub fn bit(&self, idx: usize) -> bool {
        let limb = idx / 64;
        let shift = idx % 64;
        ((self.limbs[limb] >> shift) & 1) != 0
    }

    /// Set bit at index.
    #[inline]
    pub fn set_bit(&mut self, idx: usize) {
        let limb = idx / 64;
        let shift = idx % 64;
        self.limbs[limb] |= 1u64 << shift;
    }

    /// Left shift by 1.
    #[inline]
    pub fn shl1(&mut self) {
        let mut carry = 0u64;
        for limb in &mut self.limbs {
            let next = *limb >> 63;
            *limb = (*limb << 1) | carry;
            carry = next;
        }
    }

    /// Subtraction with borrow, returns true on underflow.
    #[inline]
    pub fn overflowing_sub_assign(&mut self, rhs: &Self) -> bool {
        let mut borrow = 0u128;
        for i in 0..4 {
            let lhs = self.limbs[i] as u128;
            let rhs_i = rhs.limbs[i] as u128;
            let sub = rhs_i + borrow;
            if lhs >= sub {
                self.limbs[i] = (lhs - sub) as u64;
                borrow = 0;
            } else {
                self.limbs[i] = ((1u128 << 64) + lhs - sub) as u64;
                borrow = 1;
            }
        }
        borrow != 0
    }

    /// 128×128 → 256-bit unsigned multiply with proper carry propagation.
    /// Uses `overflowing_add` at every stage to handle inputs near u128::MAX.
    #[inline]
    pub fn mul_u128(a: u128, b: u128) -> Self {
        let a0 = a as u64 as u128;
        let a1 = (a >> 64) as u64 as u128;
        let b0 = b as u64 as u128;
        let b1 = (b >> 64) as u64 as u128;

        assemble_u128_product(a0 * b0, a0 * b1, a1 * b0, a1 * b1)
    }

    /// Lexicographic comparison of two U256 values.
    #[inline]
    pub fn cmp_words(&self, other: &Self) -> core::cmp::Ordering {
        for i in (0..4).rev() {
            match self.limbs[i].cmp(&other.limbs[i]) {
                core::cmp::Ordering::Equal => continue,
                non_eq => return non_eq,
            }
        }
        core::cmp::Ordering::Equal
    }

    /// Greater-than-or-equal comparison.
    #[inline]
    pub fn ge(&self, other: &Self) -> bool {
        self.cmp_words(other) != core::cmp::Ordering::Less
    }

    /// Divide U256 by u64, returning (quotient, remainder).
    #[inline]
    pub fn div_rem_u64(&self, divisor: u64) -> (Self, u64) {
        let mut quo = [0u64; 4];
        let mut rem = 0u64;

        for i in (0..4).rev() {
            (quo[i], rem) = div_rem_u64_step(rem, self.limbs[i], divisor);
        }

        (Self { limbs: quo }, rem)
    }

    /// Divide U256 by u128, returning (quotient, remainder).
    #[inline]
    pub fn div_rem_u128(&self, divisor: u128) -> (Self, u128) {
        const BASE: u128 = 1u128 << 64;

        let d0 = divisor as u64;
        let d1 = (divisor >> 64) as u64;
        debug_assert!(d1 != 0);

        let shift = d1.leading_zeros();
        let (v0, v1) = if shift == 0 {
            (d0, d1)
        } else {
            (d0 << shift, (d1 << shift) | (d0 >> (64 - shift)))
        };

        let mut u = [0u64; 5];
        if shift == 0 {
            u[..4].copy_from_slice(&self.limbs);
        } else {
            let mut carry = 0u64;
            for (i, limb) in self.limbs.iter().copied().enumerate() {
                u[i] = (limb << shift) | carry;
                carry = limb >> (64 - shift);
            }
            u[4] = carry;
        }

        let mut q = [0u64; 3];
        let v1_u128 = v1 as u128;
        let v0_u128 = v0 as u128;

        for j in (0..=2).rev() {
            let u2 = u[j + 2] as u128;
            let u1 = u[j + 1] as u128;
            let u0 = u[j] as u128;

            let numerator_hat = (u2 << 64) | u1;
            let mut qhat = numerator_hat / v1_u128;
            let mut rhat = numerator_hat % v1_u128;

            if qhat >= BASE {
                qhat = BASE - 1;
                rhat = numerator_hat - qhat * v1_u128;
            }

            // Knuth D3 refinement. The helper's postcondition proves the loop
            // requires at most two decrements and keeps every expression in
            // u128. Once rhat >= BASE the refinement inequality is necessarily
            // false, so evaluating rhat << 64 is neither needed nor valid.
            let (refined_qhat, _) = refine_knuth_quotient_digit(qhat, rhat, u0 as u64, v0, v1);
            qhat = refined_qhat;

            let mut carry = 0u128;
            let mut borrow = 0u128;
            for (offset, v_limb) in [v0, v1].into_iter().enumerate() {
                let product = qhat * (v_limb as u128) + carry;
                carry = product >> 64;
                let sub = (product & (BASE - 1)) + borrow;
                let cur = u[j + offset] as u128;
                if cur >= sub {
                    u[j + offset] = (cur - sub) as u64;
                    borrow = 0;
                } else {
                    u[j + offset] = (BASE + cur - sub) as u64;
                    borrow = 1;
                }
            }

            let top_sub = carry + borrow;
            let top = u[j + 2] as u128;
            let mut qhat_u64 = qhat as u64;
            if top >= top_sub {
                u[j + 2] = (top - top_sub) as u64;
            } else {
                qhat_u64 = qhat_u64.wrapping_sub(1);

                let sum0 = u[j] as u128 + v0_u128;
                u[j] = sum0 as u64;
                let carry0 = sum0 >> 64;

                let sum1 = u[j + 1] as u128 + v1_u128 + carry0;
                u[j + 1] = sum1 as u64;
                let carry1 = sum1 >> 64;

                u[j + 2] = (top + BASE - top_sub + carry1) as u64;
            }

            q[j] = qhat_u64;
        }

        let rem_norm = (u[0] as u128) | ((u[1] as u128) << 64);
        let rem = if shift == 0 {
            rem_norm
        } else {
            rem_norm >> shift
        };

        (
            Self {
                limbs: [q[0], q[1], q[2], 0],
            },
            rem,
        )
    }
}

#[cfg(kani)]
mod verification {
    use super::*;

    /// Prove the carry assembler returns the four exact base-2^64 columns for
    /// every tuple of valid 64×64 partial products. Since `mul_u128` obtains
    /// those partial products with non-overflowing u128 multiplication, this
    /// is the exact 128×128→256 multiplication postcondition.
    #[kani::proof]
    fn u256_product_carry_assembly_is_exact() {
        const BASE: u128 = 1u128 << 64;
        const MASK: u128 = BASE - 1;
        const MAX_PARTIAL: u128 = MASK * MASK;

        let p00: u128 = kani::any();
        let p01: u128 = kani::any();
        let p10: u128 = kani::any();
        let p11: u128 = kani::any();
        kani::assume(p00 <= MAX_PARTIAL);
        kani::assume(p01 <= MAX_PARTIAL);
        kani::assume(p10 <= MAX_PARTIAL);
        kani::assume(p11 <= MAX_PARTIAL);

        let product = assemble_u128_product(p00, p01, p10, p11);

        // Independent overflowing-add specification: add the unsplit middle
        // products and restore each captured 2^128 carry at the next column.
        let carry0 = p00 >> 64;
        let (middle_a, carry_a) = p01.overflowing_add(p10);
        let (middle, carry_b) = middle_a.overflowing_add(carry0);
        let carry1 = (middle >> 64) + (u128::from(carry_a) + u128::from(carry_b)) * BASE;
        let upper = p11 + carry1;

        assert_eq!(product.limbs[0], (p00 & MASK) as u64);
        assert_eq!(product.limbs[1], (middle & MASK) as u64);
        assert_eq!(product.limbs[2], (upper & MASK) as u64);
        assert_eq!(product.limbs[3], (upper >> 64) as u64);
    }

    /// Prove the concatenated dividend of one long-division step is strictly
    /// below `divisor * 2^64`, so the quotient digit fits exactly in one limb.
    /// Rust's unsigned `/` and `%` guarantees then supply the reconstruction
    /// identity and proper-remainder property for the production step.
    #[kani::proof]
    fn u256_div_rem_u64_digit_fits_one_limb() {
        let remainder: u64 = kani::any();
        let limb: u64 = kani::any();
        let divisor: u64 = kani::any();
        kani::assume(divisor != 0);
        kani::assume(remainder < divisor);

        let current = ((remainder as u128) << 64) | limb as u128;
        let exclusive_bound = (divisor as u128) << 64;
        assert!(current < exclusive_bound);
    }

    /// Prove two normalized Knuth-D3 correction transitions necessarily move
    /// `rhat` across the base. The production loop therefore performs at most
    /// two decrements; every transition is also free of underflow and overflow
    /// over its complete precondition domain.
    #[kani::proof]
    fn knuth_d3_two_corrections_cross_the_base() {
        const BASE: u128 = 1u128 << 64;

        let qhat: u64 = kani::any();
        let rhat: u128 = kani::any();
        let v1: u64 = kani::any();
        kani::assume(qhat > 0);
        kani::assume(rhat < BASE);
        kani::assume(v1 >= 1u64 << 63);

        let (qhat1, rhat1) = advance_knuth_refinement(u128::from(qhat), rhat, v1);
        assert!(rhat1 < 2 * BASE);
        assert_eq!(qhat1 + 1, u128::from(qhat));

        if rhat1 < BASE && qhat1 > 0 {
            let (qhat2, rhat2) = advance_knuth_refinement(qhat1, rhat1, v1);
            assert!(rhat2 >= BASE);
            assert!(rhat2 < 2 * BASE);
            assert_eq!(qhat2 + 2, u128::from(qhat));
        }
    }
}

// ============================================================
// Gauss-Legendre quadrature nodes and weights (test-only)
// Generated by scripts/gauss_legendre_coefficients.py — DO NOT EDIT
// ============================================================

#[cfg(test)]
pub const GL3_NODES: [i128; 3] = [-774596669241, 0, 774596669241];
#[cfg(test)]
pub const GL3_WEIGHTS: [i128; 3] = [555555555556, 888888888889, 555555555556];

#[cfg(test)]
pub const GL5_NODES: [i128; 5] = [-906179845939, -538469310106, 0, 538469310106, 906179845939];
#[cfg(test)]
pub const GL5_WEIGHTS: [i128; 5] = [
    236926885056,
    478628670499,
    568888888889,
    478628670499,
    236926885056,
];

#[cfg(test)]
pub const GL7_NODES: [i128; 7] = [
    -949107912343,
    -741531185599,
    -405845151377,
    0,
    405845151377,
    741531185599,
    949107912343,
];
#[cfg(test)]
pub const GL7_WEIGHTS: [i128; 7] = [
    129484966169,
    279705391489,
    381830050505,
    417959183673,
    381830050505,
    279705391489,
    129484966169,
];

// ============================================================
// Trigonometric functions
// Generated by scripts/trig_coefficients.py — DO NOT EDIT coefficients
// ============================================================

pub const PI_SCALE: i128 = 3_141_592_653_590;
pub const PI_OVER_2_SCALE: i128 = 1_570_796_326_795;
pub const PI_OVER_4_SCALE: i128 = 785_398_163_397;
pub const TWO_PI_SCALE: i128 = 6_283_185_307_180;
/// Sub-ULP residual of TWO_PI_SCALE: 2π·SCALE = TWO_PI_SCALE + TWO_PI_LO/SCALE.
/// 2π·1e12 = 6283185307179.5864769253…, so TWO_PI_SCALE overshoots by
/// 0.413523074713… ULP. Used by the trig range reduction (Cody-Waite).
pub const TWO_PI_LO: i128 = -413_523_074_713;

// sin(x)/x = SIN_C1 + SIN_C3·u + ... + SIN_C11·u⁵  where u = x²
pub const SIN_C1: i128 = 1_000_000_000_000;
pub const SIN_C3: i128 = -166_666_666_666;
pub const SIN_C5: i128 = 8_333_333_325;
pub const SIN_C7: i128 = -198_412_638;
pub const SIN_C9: i128 = 2_755_535;
pub const SIN_C11: i128 = -24_761;

// cos(x) = COS_C0 + COS_C2·u + ... + COS_C10·u⁵  where u = x²
pub const COS_C0: i128 = 1_000_000_000_000;
pub const COS_C2: i128 = -499_999_999_995;
pub const COS_C4: i128 = 41_666_666_560;
pub const COS_C6: i128 = -1_388_888_100;
pub const COS_C8: i128 = 24_799_027;
pub const COS_C10: i128 = -271_792;

// ============================================================
// Abramowitz-Stegun 26.2.17 constants (all scaled by 10^12)
// ============================================================

pub const INV_SQRT_2PI: i128 = 398_942_280_401;
pub const P_CDF: i128 = 231_641_900_000;
pub const B1: i128 = 319_381_530_000;
pub const B2: i128 = -356_563_782_000;
pub const B3: i128 = 1_781_477_937_000;
pub const B4: i128 = -1_821_255_978_000;
pub const B5: i128 = 1_330_274_429_000;

// ============================================================
// Piecewise Chebyshev CDF polynomial constants
// Overall max error (i128 vs mpmath): 4.2058e-11
// Both value AND derivative continuous at piece boundaries
// Generated by scripts/phi_coefficients.py — DO NOT EDIT coefficients
// ============================================================

// Interval boundaries (scaled by 1e12)
pub const POLY_I0_HI: i128 = 500_000_000_000; // 0.5
pub const POLY_I1_HI: i128 = 1_500_000_000_000; // 1.5
pub const POLY_I2_HI: i128 = 3_000_000_000_000; // 3.0
pub const POLY_I3_HI: i128 = 5_000_000_000_000; // 5.0
                                                // Piece 4 upper = 8.0 (handled by outer clamp)

// Interval midpoints and half-widths (scaled by 1e12)
pub const POLY_I0_MID: i128 = 250_000_000_000; // 0.25
pub const POLY_I0_HW: i128 = 250_000_000_000; // 0.25
pub const POLY_I1_MID: i128 = 1_000_000_000_000; // 1.0
pub const POLY_I1_HW: i128 = 500_000_000_000; // 0.5
pub const POLY_I2_MID: i128 = 2_250_000_000_000; // 2.25
pub const POLY_I2_HW: i128 = 750_000_000_000; // 0.75
pub const POLY_I3_MID: i128 = 4_000_000_000_000; // 4.0
pub const POLY_I3_HW: i128 = 1_000_000_000_000; // 1.0
pub const POLY_I4_MID: i128 = 6_500_000_000_000; // 6.5
pub const POLY_I4_HW: i128 = 1_500_000_000_000; // 1.5

// Piece 0: [0.0, 0.5] — CDF coefficients
pub const POLY_I0: [i128; 12] = [
    598_706_325_683, //  5.987063256829237e-01
    96_667_029_201,  //  9.666702920071260e-02
    -3_020_844_663,  // -3.020844662517761e-03
    -944_013_957,    // -9.440139570424519e-04
    46_217_350,      //  4.621734992857573e-05
    8_272_414,       //  8.272413997806211e-06
    -471_315,        // -4.713153491225307e-07
    -57_343,         // -5.734258314923039e-08
    3_604,           //  3.603721234054478e-09
    323,             //  3.234603705176419e-10
    -22,             // -2.170007671476414e-11
    -2,              // -1.538567210778393e-12
];

// Piece 1: [0.5, 1.5] — CDF coefficients
pub const POLY_I1: [i128; 12] = [
    841_344_746_069, //  8.413447460686025e-01
    120_985_362_260, //  1.209853622596282e-01
    -30_246_340_569, // -3.024634056900583e-02
    -2,              // -2.236691276352205e-12
    1_260_264_237,   //  1.260264237054394e-03
    -126_026_396,    // -1.260263959626473e-04
    -31_506_801,     // -3.150680066623272e-05
    6_001_162,       //  6.001162269659371e-06
    469_221,         //  4.692208922432376e-07
    -171_725,        // -1.717249631795276e-07
    -2_154,          // -2.154262323992848e-09
    3_430,           //  3.429866214864335e-09
];

// Piece 2: [1.5, 3.0] — CDF coefficients
pub const POLY_I2: [i128; 12] = [
    987_775_527_333, //  9.877755273331932e-01
    23_804_738_909,  //  2.380473890943884e-02
    -20_085_247_612, // -2.008524761161622e-02
    9_066_257_242,   //  9.066257241627936e-03
    -1_941_844_732,  // -1.941844731917149e-03
    -109_591_587,    // -1.095915869339019e-04
    176_499_579,     //  1.764995788966163e-04
    -35_211_549,     // -3.521154860402188e-05
    -3_281_426,      // -3.281425900528601e-06
    2_537_867,       //  2.537867474127547e-06
    -202_799,        // -2.027988570099364e-07
    -79_270,         // -7.927042042202421e-08
];

// Piece 3: [3.0, 5.0] — CDF coefficients
pub const POLY_I3: [i128; 12] = [
    999_968_328_754, //  9.999683287542152e-01
    133_830_262,     //  1.338302623841254e-04
    -267_660_226,    // -2.676602259835505e-04
    334_574_933,     //  3.345749333370029e-04
    -289_967_550,    // -2.899675500658869e-04
    181_788_193,     //  1.817881928678350e-04
    -82_521_635,     // -8.252163511211854e-05
    25_519_844,      //  2.551984413811646e-05
    -3_928_611,      // -3.928611463994631e-06
    -756_992,        // -7.569916104583943e-07
    656_918,         //  6.569175352291424e-07
    -150_548,        // -1.505483734178183e-07
];

// Piece 4: [5.0, 8.0] — CDF coefficients
pub const POLY_I4: [i128; 12] = [
    999_999_999_967, //  9.999999999673760e-01
    392,             //  3.915356888782659e-10
    -2_500,          // -2.499651882177424e-09
    6_511,           //  6.511131325566934e-09
    -7_968,          // -7.967599658734941e-09
    22_918,          //  2.291759574399746e-08
    -64_874,         // -6.487409521149954e-08
    57_376,          //  5.737576786733830e-08
    4_139,           //  4.138807292604606e-09
    19_426,          //  1.942643101196347e-08
    -72_091,         // -7.209080057271880e-08
    36_710,          //  3.670970644971754e-08
];

// ============================================================
// V2 HP CDF: 6 polynomial pieces + Mills ratio tail (x >= 5)
// Rounding Horner (fp_mul_hp_i already rounds), boundary-constrained,
// coordinate-descent optimized. Max 5 ULP.
// Generated by scripts/phi_hp_refit.py
// ============================================================

pub const POLY_HP_V2_I0_HI: i128 = 500_000_000_000_000;
pub const POLY_HP_V2_I1_HI: i128 = 1_500_000_000_000_000;
pub const POLY_HP_V2_I2A_HI: i128 = 2_250_000_000_000_000;
pub const POLY_HP_V2_I2B_HI: i128 = 3_000_000_000_000_000;
pub const POLY_HP_V2_I3A_HI: i128 = 4_000_000_000_000_000;

pub const POLY_HP_V2_I0_MID: i128 = 250_000_000_000_000;
pub const POLY_HP_V2_I0_HW: i128 = 250_000_000_000_000;
pub const POLY_HP_V2_I1_MID: i128 = 1_000_000_000_000_000;
pub const POLY_HP_V2_I1_HW: i128 = 500_000_000_000_000;
pub const POLY_HP_V2_I2A_MID: i128 = 1_875_000_000_000_000;
pub const POLY_HP_V2_I2A_HW: i128 = 375_000_000_000_000;
pub const POLY_HP_V2_I2B_MID: i128 = 2_625_000_000_000_000;
pub const POLY_HP_V2_I2B_HW: i128 = 375_000_000_000_000;
pub const POLY_HP_V2_I3A_MID: i128 = 3_500_000_000_000_000;
pub const POLY_HP_V2_I3A_HW: i128 = 500_000_000_000_000;
pub const POLY_HP_V2_I3B_MID: i128 = 4_500_000_000_000_000;
pub const POLY_HP_V2_I3B_HW: i128 = 500_000_000_000_000;

// I0: [0.0,0.5] deg=13 max=4 ULP
pub const POLY_HP_V2_I0: [i128; 14] = [
    598706325682924,
    96667029200713,
    -3020844662519,
    -944013957034,
    46217349936,
    8272413929,
    -471315375,
    -57342341,
    3603759,
    323049,
    -21723,
    -1210,
    4,
    -100,
];
// I1: [0.5,1.5] deg=13 max=4 ULP
pub const POLY_HP_V2_I1: [i128; 14] = [
    841344746068544,
    120985362259572,
    -30246340565023,
    -33,
    1260264191835,
    -126026418616,
    -31506612696,
    6001256127,
    468867241,
    -171906606,
    -1847116,
    3593915,
    -100208,
    -55795,
];
// I2A: [1.5,2.25] deg=15 max=4 ULP
pub const POLY_HP_V2_I2A: [i128; 16] = [
    969603638234739,
    25794853435008,
    -9068503160759,
    1520863550899,
    -54796253049,
    -24374992132,
    3883873087,
    17940493,
    -60092044,
    4454506,
    433779,
    -86366,
    3277,
    5989,
    -982,
    -1491,
];
// I2B: [2.25,3.0] deg=15 max=5 ULP
pub const POLY_HP_V2_I2B: [i128; 16] = [
    995667551636987,
    4771568098811,
    -2348506173644,
    658769960927,
    -107076056464,
    7184669478,
    828940223,
    -236847082,
    16656613,
    1411389,
    -351615,
    22524,
    5391,
    -5946,
    -829,
    1606,
];
// I3A: [3.0,4.0] deg=17 max=5 ULP
pub const POLY_HP_V2_I3A: [i128; 18] = [
    999767370920965,
    436341347522,
    -381798679096,
    204535006650,
    -73575786827,
    18081462775,
    -2821236035,
    167169387,
    39009455,
    -11645489,
    1152149,
    49566,
    -6912,
    5795,
    -13949,
    -1724,
    3518,
    416,
];
// I3B: [4.0,5.0] deg=17 max=5 ULP
pub const POLY_HP_V2_I3B: [i128; 18] = [
    999996602326874,
    7991870552,
    -8990854373,
    6410146188,
    -3231088291,
    1213608926,
    -347400516,
    75547720,
    -11941734,
    1139090,
    6868,
    -3204,
    4445,
    -25983,
    39,
    15785,
    -15,
    -3944,
];

// ---- Derivative coefficients for norm_pdf_poly ----
// d[k] = (k+1) * c[k+1], 11 coefficients per piece
// norm_pdf_poly(x) = horner_10(POLY_DERIV_Ik, t) / hw  (chain rule)

// Piece 0: [0.0, 0.5] — PDF derivative
pub const POLY_DERIV_I0: [i128; 11] = [
    96_667_029_201, //  9.666702920071260e-02
    -6_041_689_325, // -6.041689325035521e-03
    -2_832_041_871, // -2.832041871127356e-03
    184_869_400,    //  1.848693997143029e-04
    41_362_070,     //  4.136206998903105e-05
    -2_827_892,     // -2.827892094735184e-06
    -401_398,       // -4.013980820446127e-07
    28_830,         //  2.882976987243583e-08
    2_911,          //  2.911143334658777e-09
    -217,           // -2.170007671476414e-10
    -17,            // -1.692423931856233e-11
];

// Piece 1: [0.5, 1.5] — PDF derivative
pub const POLY_DERIV_I1: [i128; 11] = [
    120_985_362_260, //  1.209853622596282e-01
    -60_492_681_138, // -6.049268113801166e-02
    -7,              // -6.710073829056617e-12
    5_041_056_948,   //  5.041056948217576e-03
    -630_131_980,    // -6.301319798132367e-04
    -189_040_804,    // -1.890408039973963e-04
    42_008_136,      //  4.200813588761559e-05
    3_753_767,       //  3.753767137945901e-06
    -1_545_525,      // -1.545524668615748e-06
    -21_543,         // -2.154262323992847e-08
    37_729,          //  3.772852836350769e-08
];

// Piece 2: [1.5, 3.0] — PDF derivative
pub const POLY_DERIV_I2: [i128; 11] = [
    23_804_738_909,  //  2.380473890943884e-02
    -40_170_495_223, // -4.017049522323244e-02
    27_198_771_725,  //  2.719877172488381e-02
    -7_767_378_928,  // -7.767378927668597e-03
    -547_957_935,    // -5.479579346695095e-04
    1_058_997_473,   //  1.058997473379698e-03
    -246_480_840,    // -2.464808402281532e-04
    -26_251_407,     // -2.625140720422881e-05
    22_840_807,      //  2.284080726714792e-05
    -2_027_989,      // -2.027988570099364e-06
    -871_975,        // -8.719746246422663e-07
];

// Piece 3: [3.0, 5.0] — PDF derivative
pub const POLY_DERIV_I3: [i128; 11] = [
    133_830_262,    //  1.338302623841254e-04
    -535_320_452,   // -5.353204519671009e-04
    1_003_724_800,  //  1.003724800011009e-03
    -1_159_870_200, // -1.159870200263548e-03
    908_940_964,    //  9.089409643391750e-04
    -495_129_811,   // -4.951298106727112e-04
    178_638_909,    //  1.786389089668152e-04
    -31_428_892,    // -3.142889171195705e-05
    -6_812_924,     // -6.812924494125549e-06
    6_569_175,      //  6.569175352291424e-06
    -1_656_032,     // -1.656032107596001e-06
];

// Piece 4: [5.0, 8.0] — PDF derivative
pub const POLY_DERIV_I4: [i128; 11] = [
    392,      //  3.915356888782659e-10
    -4_999,   // -4.999303764354848e-09
    19_533,   //  1.953339397670080e-08
    -31_870,  // -3.187039863493976e-08
    114_588,  //  1.145879787199873e-07
    -389_245, // -3.892445712689973e-07
    401_630,  //  4.016303750713682e-07
    33_110,   //  3.311045834083685e-08
    174_838,  //  1.748378791076712e-07
    -720_908, // -7.209080057271879e-07
    403_807,  //  4.038067709468929e-07
];

// ============================================================
// HP (1e15) Polynomial constants for norm_cdf_poly_hp
// Generated by scripts/fit_cdf_hp.py — DO NOT EDIT
// ============================================================

// Interval boundaries (scaled by 1e15)
pub const POLY_HP_I0_HI: i128 = 500_000_000_000_000; // 0.5
pub const POLY_HP_I1_HI: i128 = 1_500_000_000_000_000; // 1.5
pub const POLY_HP_I2A_HI: i128 = 2_250_000_000_000_000; // 2.25
pub const POLY_HP_I2B_HI: i128 = 3_000_000_000_000_000; // 3.0
pub const POLY_HP_I3_HI: i128 = 5_000_000_000_000_000; // 5.0
pub const POLY_HP_I4_HI: i128 = 8_000_000_000_000_000; // 8.0

// Interval midpoints and half-widths (scaled by 1e15)
pub const POLY_HP_I0_MID: i128 = 250_000_000_000_000; // 0.25
pub const POLY_HP_I0_HW: i128 = 250_000_000_000_000; // 0.25
pub const POLY_HP_I1_MID: i128 = 1_000_000_000_000_000; // 1.0
pub const POLY_HP_I1_HW: i128 = 500_000_000_000_000; // 0.5
pub const POLY_HP_I2A_MID: i128 = 1_875_000_000_000_000; // 1.875
pub const POLY_HP_I2A_HW: i128 = 375_000_000_000_000; // 0.375
pub const POLY_HP_I2B_MID: i128 = 2_625_000_000_000_000; // 2.625
pub const POLY_HP_I2B_HW: i128 = 375_000_000_000_000; // 0.375
pub const POLY_HP_I3_MID: i128 = 4_000_000_000_000_000; // 4.0
pub const POLY_HP_I3_HW: i128 = 1_000_000_000_000_000; // 1.0
pub const POLY_HP_I4_MID: i128 = 6_500_000_000_000_000; // 6.5
pub const POLY_HP_I4_HW: i128 = 1_500_000_000_000_000; // 1.5

// HP CDF Polynomial coefficients (scaled by 1e15)
pub const POLY_HP_I0: [i128; 11] = [
    598_706_325_682_924,
    96_667_029_200_696,
    -3_020_844_662_516,
    -944_013_956_714,
    46_217_349_913,
    8_272_412_155,
    -471_315_308,
    -57_338_363,
    3_603_676,
    319_234,
    -21_682,
];

pub const POLY_HP_I1: [i128; 14] = [
    841_344_746_068_543,
    120_985_362_259_573,
    -30_246_340_564_937,
    -42,
    1_260_264_190_911,
    -126_026_418_615,
    -31_506_609_021,
    6_001_256_288,
    468_860_551,
    -171_907_068,
    -1_841_451,
    3_594_401,
    -102_019,
    -55_972,
];

pub const POLY_HP_I2A: [i128; 16] = [
    969_603_638_234_739,
    25_794_853_435_009,
    -9_068_503_160_746,
    1_520_863_550_917,
    -54_796_253_230,
    -24_374_992_417,
    3_883_874_172,
    17_942_200,
    -60_095_146,
    4_449_630,
    438_328,
    -79_214,
    -32,
    787,
    -37,
    -5,
];

pub const POLY_HP_I2B: [i128; 16] = [
    995_667_551_636_987,
    4_771_568_098_812,
    -2_348_506_173_634,
    658_769_960_908,
    -107_076_056_623,
    7_184_669_785,
    828_941_175,
    -236_848_923,
    16_653_893,
    1_416_649,
    -347_626,
    14_809,
    2_490,
    -335,
    0,
    3,
];

pub const POLY_HP_I3: [i128; 18] = [
    999_968_328_758_167,
    133_830_225_765,
    -267_660_451_532,
    334_575_564_403,
    -289_965_489_107,
    181_786_056_822,
    -82_528_639_763,
    25_518_023_995,
    -3_916_655_477,
    -740_177_127,
    644_211_338,
    -173_713_498,
    9_112_129,
    9_466_972,
    -3_317_226,
    279_600,
    129_868,
    -36_902,
];

// Dedicated HP PDF polynomial for the tail interval [5, 8].
// Used with an 8-deep Mills-ratio continued fraction to recover Φ(x) there.
pub const POLY_HP_I4_PDF: [i128; 18] = [
    266_955,
    -2_602_817,
    12_388_515,
    -38_310_277,
    86_409_890,
    -151_264_254,
    213_430_049,
    -248_623_071,
    242_831_742,
    -201_047_855,
    141_817_893,
    -84_256_548,
    41_133_871,
    -16_731_484,
    5_765_684,
    -839_426,
    -684_839,
    315_977,
];

// ============================================================
// BsIntermediates struct
// ============================================================

/// Shared intermediates for BS calculations.
/// Computes d1, d2, phi(d1), phi(d2), discount, sqrt_t, sigma_sqrt_t.
pub struct BsIntermediates {
    pub d1: i128,
    pub d2: i128,
    pub phi_d1: i128,       // Φ(d1)
    pub phi_d2: i128,       // Φ(d2)
    pub phi_neg_d1: i128,   // Φ(-d1)
    pub phi_neg_d2: i128,   // Φ(-d2)
    pub pdf_d1: i128,       // φ(d1) = norm_pdf(d1)
    pub discount: i128,     // e^(-rT)
    pub sqrt_t: i128,       // √T
    pub sigma_sqrt_t: i128, // σ√T
    pub k_disc: i128,       // K × e^(-rT)
}

// ============================================================
// BsFull struct
// ============================================================

/// Complete Black-Scholes output: prices and all Greeks at SCALE.
///
/// Returned by [`bs_full_hp`](crate::bs_full_hp) and related functions.
/// All fields are scaled by 1e12.
pub struct BsFull {
    /// Call option price (unsigned, at SCALE).
    pub call: u128,
    /// Put option price (unsigned, at SCALE).
    pub put: u128,
    /// Call delta: dC/dS (signed, at SCALE). Range [0, SCALE].
    pub call_delta: i128,
    /// Put delta: dP/dS (signed, at SCALE). Range [-SCALE, 0].
    pub put_delta: i128,
    /// Gamma: d²V/dS² (signed, at SCALE). Same for calls and puts.
    pub gamma: i128,
    /// Vega: dV/d(sigma) (signed, at SCALE). Same for calls and puts.
    pub vega: i128,
    /// Call theta: dC/dT (signed, at SCALE). Typically negative.
    pub call_theta: i128,
    /// Put theta: dP/dT (signed, at SCALE).
    pub put_theta: i128,
    /// Call rho: dC/dr (signed, at SCALE).
    pub call_rho: i128,
    /// Put rho: dP/dr (signed, at SCALE).
    pub put_rho: i128,
}

// ============================================================
// Univariate Padé [4/4] IV guess (Brenner-Subrahmanyam)
// ============================================================

// β = c·√(2π), z = x/β, σ√T ≈ β × P(z)/Q(z)
// Fitted by scripts/remez_univariate_iv.py (quick grid, 200 DE iterations)
// Max error: 0.102 σ√T on 9609-point grid
pub const SQRT_2PI_IV: i128 = 2_506_628_274_631; // √(2π) × SCALE
pub const PADE_P0: i128 = 1_030_712_890_981;
pub const PADE_P1: i128 = -3_523_716_038_286;
pub const PADE_P2: i128 = 2_286_295_838_170;
pub const PADE_P3: i128 = 148_646_789_634;
pub const PADE_P4: i128 = 72_825_165_192;
pub const PADE_Q1: i128 = -2_282_582_703_163;
pub const PADE_Q2: i128 = -135_934_203_446;
pub const PADE_Q3: i128 = -77_424_332_735;
pub const PADE_Q4: i128 = 446_572_216;

// ============================================================
// Li rational polynomial constants for implied volatility
// ============================================================

// Linear coefficients (Li 2006, original paper values at SCALE)
pub const LI_P1: i128 = -969_271_876_255; // -0.969271876255
pub const LI_P2: i128 = 97_428_338_274; //  0.097428338274
pub const LI_P3: i128 = 1_750_081_126_685; //  1.750081126685

// Numerator: n[k] for monomials x^{i_k} · (√c)^{j_k}
// k:   1  2  3  4  5  6  7  8  9 10 11 12 13 14
// i: [ 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4]
// j: [ 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0]
pub const LI_N: [i128; 14] = [
    -68_098_378_725,     // n1:  x⁰(√c)¹
    440_639_436_211,     // n2:  x¹(√c)⁰
    -263_473_754_689,    // n3:  x⁰(√c)²
    -5_792_537_721_792,  // n4:  x¹(√c)¹
    -5_267_481_008_429,  // n5:  x²(√c)⁰
    4_714_393_825_758,   // n6:  x⁰(√c)³
    3_529_944_137_559,   // n7:  x¹(√c)²
    -23_636_495_876_611, // n8:  x²(√c)¹
    -9_020_361_771_283,  // n9:  x³(√c)⁰
    14_749_084_301_452,  // n10: x⁰(√c)⁴
    -32_570_660_102_526, // n11: x¹(√c)³
    76_398_155_779_133,  // n12: x²(√c)²
    41_855_161_781_749,  // n13: x³(√c)¹
    -12_150_611_865_704, // n14: x⁴(√c)⁰
];

// Denominator: m[k] with same monomial indexing
pub const LI_M: [i128; 14] = [
    6_268_456_292_246,    // m1
    -6_284_840_445_036,   // m2
    30_068_281_276_567,   // m3
    -11_780_036_995_036,  // m4
    -2_310_966_989_723,   // m5
    -11_473_184_324_152,  // m6
    -230_101_682_610_568, // m7
    86_127_219_899_668,   // m8
    3_730_181_294_225,    // m9
    -13_954_993_561_151,  // m10
    261_950_288_864_225,  // m11
    20_090_690_444_187,   // m12
    -50_117_067_019_539,  // m13
    13_723_711_519_422,   // m14
];

// ============================================================
// AS241 inverse normal CDF coefficients (all scaled by 1e12)
// Reference: Applied Statistics algorithm AS241 (1988)
// 3-branch rational polynomial, ~12 significant digits at SCALE.
// Generated by scripts/jaeckel_coefficients.py — DO NOT EDIT
// ============================================================

// Jäckel "Let's Be Rational" solver constants
pub const SQRT_ONE_OVER_THREE: i128 = 577_350_269_190;
pub const TWO_PI_OVER_SQRT_TWENTY_SEVEN: i128 = 1_209_199_576_156;
pub const PI_OVER_SIX: i128 = 523_598_775_598;
pub const SQRT_PI_OVER_TWO: i128 = 1_253_314_137_316;
pub const TWO_PI_SCALED: i128 = 6_283_185_307_180;
pub const SQRT_THREE_SCALED: i128 = 1_732_050_807_569;

pub const AS241_SPLIT1: i128 = 425_000_000_000; // 0.425
pub const AS241_CONST1: i128 = 180_625_000_000; // 0.180625
pub const AS241_SPLIT2: i128 = 5_000_000_000_000; // 5.0
pub const AS241_CONST2: i128 = 1_600_000_000_000; // 1.6

// Branch 1 numerator: P close to 0.5
pub const AS241_A: [i128; 8] = [
    3_387_132_872_796,      // A0
    133_141_667_891_784,    // A1
    1_971_590_950_306_551,  // A2
    13_731_693_765_509_461, // A3
    45_921_953_931_549_871, // A4
    67_265_770_927_008_701, // A5
    33_430_575_583_588_128, // A6
    2_509_080_928_730_123,  // A7
];

// Branch 1 denominator
pub const AS241_B: [i128; 7] = [
    42_313_330_701_601,     // B1
    687_187_007_492_058,    // B2
    5_394_196_021_424_751,  // B3
    21_213_794_301_586_596, // B4
    39_307_895_800_092_711, // B5
    28_729_085_735_721_943, // B6
    5_226_495_278_852_855,  // B7
];

// Branch 2 numerator: P not close to 0, 0.5, or 1
pub const AS241_C: [i128; 8] = [
    1_423_437_110_750, // C0
    4_630_337_846_157, // C1
    5_769_497_221_461, // C2
    3_647_848_324_763, // C3
    1_270_458_252_452, // C4
    241_780_725_177,   // C5
    22_723_844_989,    // C6
    774_545_014,       // C7
];

// Branch 2 denominator
pub const AS241_D: [i128; 7] = [
    2_053_191_626_638, // D1
    1_676_384_830_184, // D2
    689_767_334_985,   // D3
    148_103_976_427,   // D4
    15_198_666_564,    // D5
    547_593_808,       // D6
    1_051,             // D7
];

// Branch 3 numerator: P very close to 0 or 1
pub const AS241_E: [i128; 8] = [
    6_657_904_643_501, // E0
    5_463_784_911_164, // E1
    1_784_826_539_917, // E2
    296_560_571_829,   // E3
    26_532_189_527,    // E4
    1_242_660_947,     // E5
    27_115_556,        // E6
    201_033,           // E7
];

// Branch 3 denominator
pub const AS241_F: [i128; 7] = [
    599_832_206_556, // F1
    136_929_880_923, // F2
    14_875_361_291,  // F3
    786_869_131,     // F4
    18_463_183,      // F5
    142_151,         // F6
    0,               // F7 (2.04e-15, rounds to 0)
];

// AS241 Branch 2b/2c split points and coefficients (5-branch inverse CDF).
pub const AS241_SPLIT_2B: i128 = 3_000_000_000_000; // 3.0
pub const AS241_CENTER_2B: i128 = 3_500_000_000_000; // 3.5
pub const AS241_SPLIT_2C: i128 = 4_000_000_000_000; // 4.0
pub const AS241_CENTER_2C: i128 = 4_500_000_000_000; // 4.5

pub const AS241_G: [i128; 8] = [
    4_426_662_374_924,
    88_005_498_751,
    -41_160_164_851,
    39_126_451_267,
    -58_494_244_713,
    66_910_003_962,
    33_959_185_442,
    3_528_770_658,
];
pub const AS241_H: [i128; 7] = [
    -321_373_526_953,
    105_028_553_616,
    -29_555_849_092,
    -2_051_462_415,
    15_427_578_835,
    2_488_089_624,
    293_654,
];

pub const AS241_I: [i128; 8] = [
    5_920_458_342_163,
    44_499_785_186,
    -13_437_857_419,
    10_089_893_956,
    -18_863_036_526,
    29_961_578_652,
    10_532_292_799,
    830_578_919,
];
pub const AS241_J: [i128; 7] = [
    -242_472_999_986,
    60_245_461_302,
    -14_160_150_896,
    615_344_034,
    4_831_648_864,
    586_268_402,
    34_994,
];

/// Step size for 16-entry ln table: SCALE / 16.
pub const LN_TABLE_STEP: u128 = SCALE / 16;

/// Half step for midpoint computation: SCALE / 32.
pub const LN_TABLE_HALF_STEP: u128 = SCALE / 32;

/// HP equivalents.
pub const LN_TABLE_HP_STEP: u128 = SCALE_HP_U / 16;
pub const LN_TABLE_HP_HALF_STEP: u128 = SCALE_HP_U / 32;

/// 16-entry lookup table for table-assisted ln.
/// Entry j = round(ln(1 + (2j+1)/32) × SCALE).
/// Midpoint of interval j is SCALE + (2j+1) × SCALE/32.
pub const LN_TABLE_16: [i128; 16] = [
    30771658667,  // j= 0: ln(1.03125)
    89612158690,  // j= 1: ln(1.09375)
    145182009844, // j= 2: ln(1.15625)
    197825743330, // j= 3: ln(1.21875)
    247836163905, // j= 4: ln(1.28125)
    295464212894, // j= 5: ln(1.34375)
    340926586971, // j= 6: ln(1.40625)
    384411698910, // j= 7: ln(1.46875)
    426084395311, // j= 8: ln(1.53125)
    466089729925, // j= 9: ln(1.59375)
    504556010752, // j=10: ln(1.65625)
    541597282433, // j=11: ln(1.71875)
    577315365035, // j=12: ln(1.78125)
    611801541106, // j=13: ln(1.84375)
    645137961374, // j=14: ln(1.90625)
    677398823592, // j=15: ln(1.96875)
];

/// Sub-ULP residuals for LN_TABLE_16.
/// true_ln(midpoint) × SCALE = LN_TABLE_16[j] + LN_TABLE_LO_16[j] / SCALE.
/// |LN_TABLE_LO_16[j]| < SCALE. Generated with mpmath at 60 decimal digits.
pub const LN_TABLE_LO_16: [i128; 16] = [
    -246311628972, // j= 0
    -312867380049, // j= 1
    497897281935,  // j= 2
    -80119637428,  // j= 3
    -418743219397, // j= 4
    -164123613318, // j= 5
    -406789694911, // j= 6
    332039734790,  // j= 7
    -99936875455,  // j= 8
    -400775441381, // j= 9
    395287058309,  // j=10
    -255628423458, // j=11
    -176395681888, // j=12
    -7096470110,   // j=13
    -415298334772, // j=14
    -193859190317, // j=15
];

/// HP 16-entry lookup table for table-assisted ln at SCALE_HP = 1e15.
pub const LN_TABLE_HP_16: [i128; 16] = [
    30771658666754,  // j= 0: ln(1.03125)
    89612158689687,  // j= 1: ln(1.09375)
    145182009844498, // j= 2: ln(1.15625)
    197825743329920, // j= 3: ln(1.21875)
    247836163904581, // j= 4: ln(1.28125)
    295464212893836, // j= 5: ln(1.34375)
    340926586970593, // j= 6: ln(1.40625)
    384411698910332, // j= 7: ln(1.46875)
    426084395310900, // j= 8: ln(1.53125)
    466089729924599, // j= 9: ln(1.59375)
    504556010752395, // j=10: ln(1.65625)
    541597282432744, // j=11: ln(1.71875)
    577315365034824, // j=12: ln(1.78125)
    611801541105993, // j=13: ln(1.84375)
    645137961373585, // j=14: ln(1.90625)
    677398823591806, // j=15: ln(1.96875)
];

/// Sub-ULP residuals for LN_TABLE_HP_16.
/// true_ln(midpoint) × SCALE_HP = LN_TABLE_HP_16[j] + LN_TABLE_HP_LO_16[j] / SCALE_HP.
pub const LN_TABLE_HP_LO_16: [i128; 16] = [
    -311628971792403, // j= 0
    132619951469378,  // j= 1
    -102718064936259, // j= 2
    -119637427928803, // j= 3
    256780602765747,  // j= 4
    -123613318093945, // j= 5
    210305089199780,  // j= 6
    39734790062481,   // j= 7
    63124544879595,   // j= 8
    224558619247505,  // j= 9
    287058308531738,  // j=10
    371576542303900,  // j=11
    -395681887938481, // j=12
    -96470110233571,  // j=13
    -298334771503865, // j=14
    140809682609997,  // j=15
];

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

    #[test]
    fn test_ln2_lo_sub_ulp() {
        assert!(
            (LN2_LO as u128) < SCALE || ((-LN2_LO) as u128) < SCALE,
            "LN2_LO={} not sub-ULP",
            LN2_LO
        );
    }

    #[test]
    fn test_ln2_hp_lo_sub_ulp() {
        assert!(
            (LN2_HP_LO.unsigned_abs()) < SCALE_HP_U,
            "LN2_HP_LO={} not sub-ULP",
            LN2_HP_LO
        );
    }

    #[test]
    fn test_ln2_lo_no_overflow() {
        for k in -127i128..=127 {
            assert!(k.checked_mul(LN2_LO).is_some(), "overflow at k={}", k);
        }
    }

    #[test]
    fn test_ln2_hp_lo_no_overflow() {
        for k in -127i128..=127 {
            assert!(k.checked_mul(LN2_HP_LO).is_some(), "overflow at k={}", k);
        }
    }

    #[test]
    fn test_ln2_correction_bounded() {
        for k in -88i128..=88 {
            let raw = k * LN2_LO;
            let correction = if raw >= 0 {
                (raw + SCALE_I / 2) / SCALE_I
            } else {
                (raw - SCALE_I / 2) / SCALE_I
            };
            assert!(
                correction.abs() <= 10,
                "Correction too large at k={}: {}",
                k,
                correction
            );
        }
    }

    #[test]
    fn test_ln_remez_coeffs_matches_individual() {
        assert_eq!(LN_REMEZ_COEFFS[0], LN_REMEZ_W0);
        assert_eq!(LN_REMEZ_COEFFS[1], LN_REMEZ_W1);
        assert_eq!(LN_REMEZ_COEFFS[7], LN_REMEZ_W7);
    }

    #[test]
    fn test_ln_remez_hp_coeffs_matches_individual() {
        assert_eq!(LN_REMEZ_HP_COEFFS[0], LN_REMEZ_HP0);
        assert_eq!(LN_REMEZ_HP_COEFFS[9], LN_REMEZ_HP9);
    }

    #[test]
    fn test_ln_table_lo_sub_ulp() {
        for j in 0..16 {
            assert!(
                LN_TABLE_LO_16[j].unsigned_abs() < SCALE,
                "LN_TABLE_LO_16[{}]={} not sub-ULP",
                j,
                LN_TABLE_LO_16[j]
            );
        }
    }

    #[test]
    fn test_ln_table_hp_lo_sub_ulp() {
        for j in 0..16 {
            assert!(
                LN_TABLE_HP_LO_16[j].unsigned_abs() < SCALE_HP_U,
                "LN_TABLE_HP_LO_16[{}]={} not sub-ULP",
                j,
                LN_TABLE_HP_LO_16[j]
            );
        }
    }

    #[test]
    fn test_ln_table_monotone() {
        for j in 1..16 {
            assert!(
                LN_TABLE_16[j] > LN_TABLE_16[j - 1],
                "LN_TABLE_16 not monotone at j={}",
                j
            );
            assert!(
                LN_TABLE_HP_16[j] > LN_TABLE_HP_16[j - 1],
                "LN_TABLE_HP_16 not monotone at j={}",
                j
            );
        }
    }
}