ported/chi_f_t_family

fn truncate_positive_integer(value: f64) -> Result<f64, WorksheetErrorCode> {
    let truncated = value.trunc();
    // Excel docs: #NUM! when deg_freedom > 1e10 (strictly greater); 1e10 itself is admitted.
    if !truncated.is_finite() || truncated < 1.0 || truncated > 1e10 {
        Err(WorksheetErrorCode::Num)
    } else {
        Ok(truncated)
    }
}

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

fn validate_nonnegative_x(x: f64) -> Result<f64, WorksheetErrorCode> {
    if !x.is_finite() || x < 0.0 {
        Err(WorksheetErrorCode::Num)
    } else {
        Ok(x)
    }
}

pub fn chisq_dist_rt_kernel(x: f64, deg_freedom: f64) -> Result<f64, WorksheetErrorCode> {
    let x = validate_nonnegative_x(x)?;
    let k = truncate_positive_integer(deg_freedom)?;
    // Same capture: CHIDIST(x,1) == CHISQ.DIST.RT(x,1) == ERFC.PRECISE(SQRT(x/2))
    // bit-exactly. Routing through the published ERFC kernel is the identified
    // graph, not a numeric fit.
    if k == 1.0 {
        erfc_of_sqrt_half_x(x)
    } else if k == 2.0 {
        // Inverse-problem identity, live Excel 16.0 b20228: CHIDIST(x,2)
        // == EXP(-x/2) == EXP(-(x/2)) on 68/68 nonnegative rows. The
        // worksheet EXP is the signed-off x87 elementary. The complementary
        // CDF 1-EXP(-x/2) is NOT this graph (9/11 on a follow-up bank).
        Ok(crate::ported::excel_numeric::excel_exp(-(x / 2.0)))
    } else {
        Ok(regularized_gamma_q(k / 2.0, x / 2.0))
    }
}

fn search_upper_bound<F>(target: f64, initial_hi: f64, f: F) -> f64
where
    F: Fn(f64) -> f64,
{
    let mut hi = initial_hi.max(1.0);
    for _ in 0..200 {
        if f(hi) >= target {
            return hi;
        }
        hi *= 2.0;
    }
    hi
}

pub fn chisq_inv_rt_kernel(probability: f64, deg_freedom: f64) -> Result<f64, WorksheetErrorCode> {
    let p = validate_probability_open_unit(probability)?;
    let k = truncate_positive_integer(deg_freedom)?;
    // Excel inverts its published right-tail surface (CHIDIST = Q) directly,
    // not P at 1-p: the complement staging carries a systematic -5..-33-ULP
    // bias on the b14 CHIINV corpus (W109). Q is decreasing in x, so invert
    // the negated forward at -p to keep the increasing-forward convention.
    let hi = search_upper_bound(-p, k, |x| -regularized_gamma_q(k / 2.0, x / 2.0));
    Ok(bisect_inverse(-p, 0.0, hi, |x| {
        -regularized_gamma_q(k / 2.0, x / 2.0)
    }))
}

pub fn f_dist_rt_kernel(x: f64, deg1: f64, deg2: f64) -> Result<f64, WorksheetErrorCode> {
    let x = validate_nonnegative_x(x)?;
    let d1 = truncate_positive_integer(deg1)?;
    let d2 = truncate_positive_integer(deg2)?;
    // W109 identified staging: x = d2/den, y = d1*F/den (accurate complement).
    let num = d1 * x;
    let den = d2 + num;
    Ok(bratio(d2 / 2.0, d1 / 2.0, d2 / den, num / den).0)
}

pub fn f_inv_rt_kernel(probability: f64, deg1: f64, deg2: f64) -> Result<f64, WorksheetErrorCode> {
    let p = validate_probability_open_unit(probability)?;
    let d1 = truncate_positive_integer(deg1)?;
    let d2 = truncate_positive_integer(deg2)?;
    // Invert the published right-tail surface (f_dist_rt's accurate complement
    // form) directly at p, not the CDF at 1-p — same principle as CHIINV
    // (W109 b14+b19: the complement staging carries a systematic small-p bias).
    // The surface is decreasing in x, so invert the negated forward at -p.
    let f = move |x: f64| {
        let num = d1 * x;
        let den = d2 + num;
        -bratio(d2 / 2.0, d1 / 2.0, d2 / den, num / den).0
    };
    let hi = search_upper_bound(-p, 1.0, f);
    Ok(bisect_inverse(-p, 0.0, hi, f))
}

fn t_cdf(x: f64, deg_freedom: f64) -> Result<f64, WorksheetErrorCode> {
    let v = truncate_positive_integer(deg_freedom)?;
    // W109 identified staging: x = df/den, y = t^2/den (accurate complement).
    let t2 = x * x;
    let den = v + t2;
    let ib = bratio(v / 2.0, 0.5, v / den, t2 / den).0;
    Ok(if x >= 0.0 { 1.0 - 0.5 * ib } else { 0.5 * ib })
}

pub fn t_dist_rt_kernel(x: f64, deg_freedom: f64) -> Result<f64, WorksheetErrorCode> {
    let x = validate_nonnegative_x(x)?;
    let v = truncate_positive_integer(deg_freedom)?;
    let t2 = x * x;
    let den = v + t2;
    Ok(0.5 * bratio(v / 2.0, 0.5, v / den, t2 / den).0)
}

pub fn t_dist_2t_kernel(x: f64, deg_freedom: f64) -> Result<f64, WorksheetErrorCode> {
    let x = validate_nonnegative_x(x)?;
    let v = truncate_positive_integer(deg_freedom)?;
    let t2 = x * x;
    let den = v + t2;
    Ok(bratio(v / 2.0, 0.5, v / den, t2 / den).0)
}

pub fn t_inv_2t_kernel(probability: f64, deg_freedom: f64) -> Result<f64, WorksheetErrorCode> {
    let p = validate_probability_open_unit(probability)?;
    let v = truncate_positive_integer(deg_freedom)?;
    if v == 1.0 {
        // Live Excel 16.0 b20326: T.INV.2T(p,1)=1/TAN(PI()*p/2) 9/9.
        // Cauchy closed form; the BETA.INV sqrt graph misses 3/9 here.
        // p=0.5 is the exact quartile 1; worksheet TAN(PI()/4) is 1, but
        // excel_tan(π/4) is 1 ULP high so 1/tan would publish 1+1 ULP.
        if p == 0.5 {
            return Ok(1.0);
        }
        let arg = std::f64::consts::PI * p / 2.0;
        return Ok(1.0 / crate::ported::tan::tan_kernel(arg));
    }
    // df>=2: Excel is SQRT(df*(1/BETA.INV(p,df/2,0.5)-1)) 135/135 vs the
    // worksheet BETA.INV, but production BETA.INV is still a few ULP off
    // Excel, so the composition is not landable yet. Keep the published
    // two-tail bratio inverter (W109 b19).
    let f = move |x: f64| {
        let t2 = x * x;
        let den = v + t2;
        -bratio(v / 2.0, 0.5, v / den, t2 / den).0
    };
    let hi = search_upper_bound(-p, 1.0, f);
    Ok(bisect_inverse(-p, 0.0, hi, f))
}