ported/normal_log_family

fn flush_subnormal(value: f64) -> f64 {
    if value.abs() < f64::MIN_POSITIVE {
        0.0
    } else {
        value
    }
}

/// G-F3: `NORMSDIST` / `NORM.S.DIST` CDF wrapper. W111 mathematical
/// discriminators require `z = RN53(RN64(|x|*RN(1/√2)))`; native binary64
/// multiplication fails 133/512 fresh public-ERFC composition controls.
/// `Q` is the published ERFC body at `z` (still an open kernel);
/// `x < 0 → RN53(0.5*Q)`, `x ≥ 0 → RN53(1-0.5*Q)`; PHI-class flush after
/// the 0.5 multiply (`NORMSDIST(-37.52) = +0` with ERFC still finite).
pub fn identified_std_normal_cdf(x: f64) -> f64 {
    if x == 0.0 {
        return 0.5;
    }
    let z = crate::ported::excel_numeric::excel_x87_mul(x.abs(), f64::from_bits(FRAC_1_SQRT_2_BITS));
    let q = match erfc_precise_kernel(z) {
        Ok(value) => value,
        Err(_) => return f64::NAN,
    };
    let half_q = 0.5 * q;
    let ns = if x.is_sign_negative() {
        half_q
    } else {
        1.0 - half_q
    };
    flush_subnormal(ns)
}

fn norm_cdf(x: f64) -> f64 {
    identified_std_normal_cdf(x)
}

fn validate_positive_sigma(sigma: f64) -> Result<(), WorksheetErrorCode> {
    if sigma <= 0.0 || !sigma.is_finite() {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(())
}

fn validate_probability_open_unit(p: f64) -> Result<(), WorksheetErrorCode> {
    if !p.is_finite() || p <= 0.0 || p >= 1.0 {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(())
}

fn inverse_standard_normal_acklam_refined(p: f64) -> f64 {
    const A: [f64; 6] = [
        -3.969_683_028_665_376e1,
        2.209_460_984_245_205e2,
        -2.759_285_104_469_687e2,
        1.383_577_518_672_69e2,
        -3.066_479_806_614_716e1,
        2.506_628_277_459_239,
    ];
    const B: [f64; 5] = [
        -5.447_609_879_822_406e1,
        1.615_858_368_580_409e2,
        -1.556_989_798_598_866e2,
        6.680_131_188_771_972e1,
        -1.328_068_155_288_572e1,
    ];
    const C: [f64; 6] = [
        -7.784_894_002_430_293e-3,
        -3.223_964_580_411_365e-1,
        -2.400_758_277_161_838,
        -2.549_732_539_343_734,
        4.374_664_141_464_968,
        2.938_163_982_698_783,
    ];
    const D: [f64; 4] = [
        7.784_695_709_041_462e-3,
        3.224_671_290_700_398e-1,
        2.445_134_137_142_996,
        3.754_408_661_907_416,
    ];
    const P_LOW: f64 = 0.02425;
    const P_HIGH: f64 = 1.0 - P_LOW;

    let x = if p < P_LOW {
        let q = (-2.0 * crt_log(p)).sqrt();
        (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
    } else if p <= P_HIGH {
        let q = p - 0.5;
        let r = q * q;
        (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q
            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
    } else {
        let q = (-2.0 * crt_log(1.0 - p)).sqrt();
        -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
    };

    // Form the lower tail directly. Subtracting erf from 1 loses it entirely
    // for small probabilities; at positive x compute the complementary error
    // before subtracting, so the Newton residual does not round through 1.
    let tail = 0.5 * libm::erfc(x.abs() / std::f64::consts::SQRT_2);
    let err = if x < 0.0 { tail - p } else { (1.0 - p) - tail };
    let density = phi_kernel(x);
    if density == 0.0 { x } else { x - err / density }
}

fn inverse_standard_normal(p: f64) -> Result<f64, WorksheetErrorCode> {
    validate_probability_open_unit(p)?;

    Ok(inverse_standard_normal_acklam_refined(p))
}

pub fn norm_s_dist_kernel(z: f64, cumulative: bool) -> Result<f64, WorksheetErrorCode> {
    if cumulative {
        Ok(norm_cdf(z))
    } else {
        crate::ported::phi_fn::phi_kernel(z)
    }
}

pub fn norm_dist_kernel(
    x: f64,
    mean: f64,
    sigma: f64,
    cumulative: bool,
) -> Result<f64, WorksheetErrorCode> {
    validate_positive_sigma(sigma)?;
    if cumulative {
        let z = (x - mean) / sigma;
        Ok(norm_cdf(z))
    } else {
        // W111 density observations distinguish all these publication points:
        // subtract and divide at RN64 then store binary64; divide the unflushed
        // exponential by sigma before multiplying the normal constant. Calling
        // worksheet PHI first loses both scale ordering and tiny intermediates.
        use crate::ported::excel_numeric::{excel_exp, excel_x87_div, excel_x87_mul, excel_x87_sub};
        let z = excel_x87_div(excel_x87_sub(x, mean), sigma);
        let square = excel_x87_mul(z, z);
        if !square.is_finite() {
            return Err(WorksheetErrorCode::Num);
        }
        let exponential = excel_exp(-(square / 2.0));
        let value = excel_x87_mul(
            excel_x87_div(exponential, sigma),
            f64::from_bits(0x3fd9884533d43651),
        );
        Ok(flush_subnormal(value))
    }
}

pub fn norm_s_inv_kernel(p: f64) -> Result<f64, WorksheetErrorCode> {
    inverse_standard_normal(p)
}

pub fn norm_inv_kernel(p: f64, mean: f64, sigma: f64) -> Result<f64, WorksheetErrorCode> {
    validate_positive_sigma(sigma)?;
    Ok(mean + sigma * inverse_standard_normal(p)?)
}

pub fn lognorm_dist_kernel(
    x: f64,
    mean: f64,
    sigma: f64,
    cumulative: bool,
) -> Result<f64, WorksheetErrorCode> {
    validate_positive_sigma(sigma)?;
    if x <= 0.0 || !x.is_finite() {
        return Err(WorksheetErrorCode::Num);
    }
    // Inverse-problem identity, live Excel 16.0 b20228: LOGNORM.DIST(x,μ,σ,TRUE)
    // == NORM.S.DIST((LN(x)-μ)/σ,TRUE) == NORMSDIST((LN(x)-μ)/σ)
    // == NORM.DIST(LN(x),μ,σ,TRUE) on 45/45 mixed (x,μ,σ) rows.
    // The full recorded banks distinguish the internal CRT logarithm at oc3066;
    // the public LN chain differs there. The earlier 45-row identity still holds.
    let z = (crt_log(x) - mean) / sigma;
    if cumulative {
        Ok(norm_cdf(z))
    } else {
        Ok(phi_kernel(z) / (x * sigma))
    }
}

pub fn lognorm_inv_kernel(p: f64, mean: f64, sigma: f64) -> Result<f64, WorksheetErrorCode> {
    validate_positive_sigma(sigma)?;
    // Same capture: LOGNORM.INV(0.5,μ,σ) == EXP(NORM.INV(0.5,μ,σ)) on 45/45.
    Ok(crate::ported::excel_numeric::excel_exp(
        mean + sigma * inverse_standard_normal(p)?,
    ))
}

pub fn confidence_norm_kernel(
    alpha: f64,
    standard_dev: f64,
    size: f64,
) -> Result<f64, WorksheetErrorCode> {
    validate_probability_open_unit(alpha)?;
    if standard_dev <= 0.0 || !standard_dev.is_finite() || !size.is_finite() {
        return Err(WorksheetErrorCode::Num);
    }
    // Excel docs: size is truncated toward integer; size < 1 -> #NUM!.
    let size_int = size.trunc();
    if size_int < 1.0 {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(inverse_standard_normal(1.0 - alpha / 2.0)? * standard_dev / size_int.sqrt())
}