ported/normal_dist_common

/// Worksheet `PHI` numerical graph observed on the reference Excel host.
/// W111 retains a fresh 9,130-row exact numerical replay; surface overflow
/// admission and contextual coercion remain separate obligations.
///
/// ```text
/// sq  = RN53(RN64(x·x))            (x87 double-rounded square)
/// t   = -(sq/2)                    (exact halving + negation)
/// e   = EXP(t)                     (the x87 exponential, `excel_exp`)
/// v   = RN53(RN64(e · RN(1/sqrt(2π))))
/// PHI = +0.0 if |v| < DBL_MIN      (live-pinned publication flush;
///        normal values publish unchanged — boundary probed at bit level)
/// ```
pub fn phi_kernel(x: f64) -> f64 {
    use crate::ported::excel_numeric::{excel_exp, excel_x87_mul};
    let sq = excel_x87_mul(x, x);
    let e = excel_exp(-(sq / 2.0));
    let v = excel_x87_mul(e, INV_SQRT_2PI);
    if v.abs() < f64::MIN_POSITIVE { 0.0 } else { v }
}

/// One Horner evaluation with every product and sum stored to binary64, as `erf_approx` needs.
fn stored_horner(coefficients: &[u64], t: f64) -> f64 {
    use crate::ported::excel_numeric::{excel_x87_add, excel_x87_mul};
    let mut value = 0.0;
    for (i, &bits) in coefficients.iter().enumerate().rev() {
        let c = f64::from_bits(bits);
        value = if i + 1 == coefficients.len() { c } else { excel_x87_add(excel_x87_mul(value, t), c) };
    }
    value
}

/// The shared ERF body (ERF, ERF.PRECISE, the below-half ERFC complement and GAUSS's tiny
/// route).
///
/// For `|x| < 0.5` Excel evaluates a (4, 5) rational `x * N(t) / D(t)` as compiled double
/// arithmetic on the x87 at precision 64: every operation is `RN53(RN64(..))`, that is, its
/// result is stored to binary64 before the next one reads it, in the order
///
/// ```text
/// t = x*x;  N = Horner(t);  D = Horner(t);  p = x*N;  y = p/D
/// ```
///
/// Excel shows this directly. Stored `N` and `D` step at edges found by bisection to one input
/// ulp, and the two `D` edges measured agree to 1e-15 only with the 64-bit register's half unit on
/// each. Of the store placements raced on 2,000 consecutive answers at each of about 20 sites up
/// to 0.5, all stores leaves 0.5% of the rows unexplained (the coefficients' last bits) and every
/// other placement more than twice that. The measured coefficients give `N1/N0 - D1/D0 = -1/3`
/// to 2e-16 and reproduce erf to a flat offset of 9e-17 up to 0.5, past which the error climbs at
/// once, which places the approximation's design interval.
///
/// Above 0.5 the magnitude is the stored-PC64 complement of the upper ERFC body. Independent
/// public calls agree on 3,328 adjacent inputs; exponent-edge controls distinguish the final
/// PC64 store from strict binary64 subtraction. Up to 4 ERFC's `F` remains a stand-in, so this measured
/// composition does not establish its last bits there. Exactly 0.5 keeps its existing boundary body.
/// Subnormal arguments cannot be written to a cell (Excel reads them back as zero).
pub fn erf_approx(x: f64) -> f64 {
    if x.is_normal() && x.abs() < 0.5 {
        use crate::ported::excel_numeric::{excel_x87_div, excel_x87_mul};
        let t = excel_x87_mul(x, x);
        let numerator = stored_horner(&ERF_N, t);
        let denominator = stored_horner(&ERF_D, t);
        return excel_x87_div(excel_x87_mul(x, numerator), denominator);
    }
    if x.abs() > 0.5 {
        use crate::ported::excel_numeric::excel_x87_sub;
        return excel_x87_sub(1.0, erfc_upper(x.abs())).copysign(x);
    }
    libm::erf(x)
}

/// The upper branch of Excel's erfc for `x >= 0.5`: `Q = E * F`.
///
/// `E` is measured. `z = RN53(RN64(x*x))` and `E = exp(-z)` are stored to binary64 before use: over
/// 220,000 Excel answers up to `x = 26.5` the residual of `Q` against this `E` times the true
/// `exp(x*x)*erfc(x)` stays within about 10 units of 2^-53, where the exact `exp(-x*x)` spreads by
/// hundreds of units. `z`'s double rounding shows directly: plain binary64 `x*x` leaves outliers of
/// up to 1,024 units.
///
/// `F` has two branches. `x = 4` itself falls below: bisected on consecutive doubles, Excel leaves the table below at the next double up.
///
/// Above 4, `F = N(x) / D(x)`, degrees (8, 9), Horner in `x`, each product, sum and the final
/// quotient stored to binary64 from precision 64, and `Q` is the stored product `E * F`. The
/// coefficients were read from Excel one stored stage at a time, from the top of each polynomial:
/// the rounding of a stored stage leaves a sawtooth in the answers, so runs of consecutive doubles
/// at geometrically spaced `x` give each stage's phase, and the next stage's phase over a long
/// baseline gives the coefficient above it. The answers' last bits then fix the low bits, band by
/// band downward from `x = 26.5`. `D(0)` reads as 1 and `N/D` follows `exp(x*x)*erfc(x)`'s
/// expansion at infinity: the top two coefficients agree between `N` and `D` to 1e-14, the next
/// differs by 1/2.
///
/// Up to 4 `F` is a stand-in with the same store discipline, and last bits differ. The stage tones
/// of the table above stop between 3.95 and 4.05. Excel's own table there is not identified.
pub fn erfc_upper(x: f64) -> f64 {
    use crate::ported::excel_numeric::{excel_exp, excel_x87_add, excel_x87_div, excel_x87_mul};
    let step = |value: f64, c: f64, by: f64| excel_x87_mul(excel_x87_add(value, c), by);
    let z = excel_x87_mul(x, x);
    let e = excel_exp(-z);
    if e == 0.0 {
        // Past about 27 the exponential is gone and the polynomials below can overflow.
        return 0.0;
    }
    let f = if x > 4.0 {
        excel_x87_div(stored_horner(&ERFC_TAIL_N, x), stored_horner(&ERFC_TAIL_D, x))
    } else {
        let mut num = excel_x87_mul(ERFC_C[8], x);
        let mut den = x;
        for i in 0..7 {
            num = step(num, ERFC_C[i], x);
            den = step(den, ERFC_D[i], x);
        }
        excel_x87_div(excel_x87_add(num, ERFC_C[7]), excel_x87_add(den, ERFC_D[7]))
    };
    excel_x87_mul(e, f)
}