ported/special_dist_family

pub fn erf_precise_kernel(x: f64) -> Result<f64, WorksheetErrorCode> {
    if !x.is_finite() {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(erf_approx(x))
}

pub fn erfc_of_sqrt_half_x(x: f64) -> Result<f64, WorksheetErrorCode> {
    erfc_precise_kernel((x / 2.0).sqrt())
}

pub fn erfc_kernel(x: f64) -> Result<f64, WorksheetErrorCode> {
    if !x.is_finite() {
        return Err(WorksheetErrorCode::Num);
    }
    // Recorded ERFC/ERFC.PRECISE tail readings publish +0 below the normal range.
    let value = if x > 0.0 && x <= 0.5 {
        1.0 - erf_approx(x)
    } else if x > 0.5 {
        crate::ported::normal_dist_common::erfc_upper(x)
    } else if x < -0.5 {
        // R30 public complement pairs. The exponent controls distinguish the PC64
        // subtraction from strict binary64. The positive upper body still owns F.
        crate::ported::excel_numeric::excel_x87_sub(
            2.0, crate::ported::normal_dist_common::erfc_upper(-x),
        )
    } else {
        libm::erfc(x)
    };
    Ok(if value.is_subnormal() { 0.0 } else { value })
}

pub fn erfc_precise_kernel(x: f64) -> Result<f64, WorksheetErrorCode> {
    erfc_kernel(x)
}

pub fn gammaln_kernel(x: f64) -> Result<f64, WorksheetErrorCode> {
    // Published GAMMALN / GAMMALN.PRECISE surface. The full positive numeric
    // graph remains under W109 identification; this admission boundary is
    // independently pinned. Current-reference Excel rejects every positive
    // binary64 subnormal with #NUM! while admitting `f64::MIN_POSITIVE`; the
    // same guard also retains the non-positive domain. GAMMA and the shared
    // internal lgamma are unaffected.
    if !x.is_finite() || x < f64::MIN_POSITIVE {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(crate::ported::excel_numeric::gammaln_excel(x))
}

pub fn weibull_dist_kernel(
    x: f64,
    alpha: f64,
    beta: f64,
    cumulative: bool,
) -> Result<f64, WorksheetErrorCode> {
    if !x.is_finite()
        || !alpha.is_finite()
        || !beta.is_finite()
        || x < 0.0
        || alpha <= 0.0
        || beta <= 0.0
    {
        return Err(WorksheetErrorCode::Num);
    }

    if x == 0.0 {
        if cumulative {
            return Ok(0.0);
        }
        return Ok(0.0);
    }

    // W109 lane-1 identification (b24 + b27/b27b, all blocks bit-exact): the
    // WEIBULL body is a legacy x87 compilation unit — every op double-rounded
    // through a spilled double local, every pow the raw chain (no shortcuts).
    //   r   = RN53(RN64(x/β))                                   (b27b D2)
    //   t   = exp(RN53(RN64(α·ln r)))                           (b27D 113/113)
    //   cdf = −expm1(−t)                                        (Kahan expm1)
    //   pdf = DR(DR(DR(α/β^α)·x^(α−1))·exp(−t))                 (b27 T3|SS,
    //         division-first association: `alpha / pow(beta, alpha) *
    //         pow(x, alpha-1) * exp(-pow(x/beta, alpha))` left-to-right)
    use crate::ported::excel_numeric::{
        excel_exp, excel_expm1_internal, excel_pow_chain, excel_x87_div, excel_x87_mul,
    };
    let ratio = excel_x87_div(x, beta);
    let power = excel_pow_chain(ratio, alpha);
    let value = if cumulative {
        -excel_expm1_internal(-power)
    } else {
        let pba = excel_pow_chain(beta, alpha);
        let px = excel_pow_chain(x, alpha - 1.0);
        let e = excel_exp(-power);
        excel_x87_mul(excel_x87_mul(excel_x87_div(alpha, pba), px), e)
    };

    if !value.is_finite() {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(value)
}