MULTINOMIAL

40 multinomial "MULTINOMIAL" (1, 255) [R "number1", R "number2"] |s, a, sv| {
    // Admission runs across every value before the numeric domain is checked.
    // Held blank cells are skipped; logicals are rejected in every origin.
    let mut numbers = Vec::new();
    values(s, a, 0, a.len(), sv, &mut |held, value| {
        if held && matches!(value.core_ref(), CoreRef::Empty) {
            return Ok(());
        }
        reject_logical(value)?;
        numbers.push(s.number_argument(0, value)?);
        Ok(())
    })?;
    let mut items = Vec::with_capacity(numbers.len());
    for n in numbers {
        if !n.is_finite() || n < 0.0 {
            return fail(WorksheetErrorCode::Num);
        }
        items.push(n.trunc() as i64);
    }
    out(s, k::multinomial::multinomial_kernel(&items))
}

pub fn factorial_of_int(n: i64) -> f64 {
    if n <= 1 {
        return 1.0;
    }
    (2..=n).fold(1.0, |acc, item| acc * item as f64)
}

fn next_down(x: f64) -> f64 {
    if x.is_nan() || x == f64::NEG_INFINITY {
        x
    } else if x == 0.0 {
        -f64::from_bits(1)
    } else if x > 0.0 {
        f64::from_bits(x.to_bits() - 1)
    } else {
        f64::from_bits(x.to_bits() + 1)
    }
}

fn step_down(mut x: f64, steps: u64) -> f64 {
    for _ in 0..steps {
        x = next_down(x);
    }
    x
}

fn ulp_distance(lhs: f64, rhs: f64) -> u64 {
    lhs.to_bits().abs_diff(rhs.to_bits())
}

fn multinomial_positive_items(items: &[i64]) -> Vec<i64> {
    items.iter().copied().filter(|item| *item > 0).collect()
}

fn positive_items_are_consecutive(positive_items: &[i64]) -> bool {
    if positive_items.len() < 2 {
        return false;
    }
    let mut sorted = positive_items.to_vec();
    sorted.sort_unstable();
    sorted.windows(2).all(|window| window[1] == window[0] + 1)
}

fn multinomial_exact_candidate(positive_items: &[i64]) -> f64 {
    let total = positive_items.iter().sum::<i64>();
    let denominator = positive_items
        .iter()
        .fold(1.0, |acc, item| acc * factorial_of_int(*item));
    factorial_of_int(total) / denominator
}

fn multinomial_log_factorial_candidate(positive_items: &[i64]) -> f64 {
    let total = positive_items.iter().sum::<i64>();
    let numerator_ln = crt_log(factorial_of_int(total));
    let denominator_ln = positive_items
        .iter()
        .map(|item| crt_log(factorial_of_int(*item)))
        .sum::<f64>();
    crt_exp(numerator_ln - denominator_ln)
}

fn lanczos_ln_gamma(x: f64) -> f64 {
    let z = x - 1.0;
    let mut acc = MULTINOMIAL_LANCZOS_COEFFS[0];
    for (index, coeff) in MULTINOMIAL_LANCZOS_COEFFS.iter().enumerate().skip(1) {
        acc += coeff / (z + index as f64);
    }

    let t = z + MULTINOMIAL_LANCZOS_G + 0.5;
    0.5 * crt_log(2.0 * std::f64::consts::PI) + (z + 0.5) * crt_log(t) - t + crt_log(acc)
}

fn multinomial_descending_log_numerator_lanczos_denominator_candidate(
    positive_items: &[i64],
) -> f64 {
    let total = positive_items.iter().sum::<i64>();
    let numerator_ln = (2..=total)
        .rev()
        .map(|value| crt_log(value as f64))
        .sum::<f64>();
    let denominator_ln = positive_items
        .iter()
        .map(|item| lanczos_ln_gamma((*item + 1) as f64))
        .sum::<f64>();
    crt_exp(numerator_ln - denominator_ln)
}

fn multinomial_lanczos_candidate(positive_items: &[i64]) -> f64 {
    let total = positive_items.iter().sum::<i64>();
    let numerator_ln = lanczos_ln_gamma((total + 1) as f64);
    let denominator_ln = positive_items
        .iter()
        .map(|item| lanczos_ln_gamma((*item + 1) as f64))
        .sum::<f64>();
    crt_exp(numerator_ln - denominator_ln)
}

fn multinomial_publication_candidate(positive_items: &[i64], exact_candidate: f64) -> f64 {
    let total = positive_items.iter().sum::<i64>();

    match positive_items.len() {
        0 | 1 => exact_candidate,
        2 => {
            if total < 5 {
                return exact_candidate;
            }

            let candidate = multinomial_log_factorial_candidate(positive_items);
            if candidate.is_finite() && candidate < exact_candidate {
                let floor = step_down(
                    exact_candidate,
                    MULTINOMIAL_TWO_POSITIVE_PUBLICATION_ULP_CLAMP,
                );
                candidate.max(floor)
            } else {
                exact_candidate
            }
        }
        3 => {
            if total < 9 || !positive_items_are_consecutive(positive_items) {
                return exact_candidate;
            }

            let candidate = multinomial_log_factorial_candidate(positive_items);
            if candidate.is_finite() && candidate < exact_candidate {
                candidate
            } else {
                exact_candidate
            }
        }
        _ => {
            if total < 10 || !positive_items_are_consecutive(positive_items) {
                return exact_candidate;
            }

            let mut best = exact_candidate;
            let mut best_distance = MULTINOMIAL_MULTI_POSITIVE_PUBLICATION_ULP_RADIUS + 1;
            for candidate in [
                multinomial_descending_log_numerator_lanczos_denominator_candidate(positive_items),
                multinomial_lanczos_candidate(positive_items),
            ] {
                if candidate.is_finite() {
                    let distance = ulp_distance(candidate, exact_candidate);
                    if distance <= MULTINOMIAL_MULTI_POSITIVE_PUBLICATION_ULP_RADIUS
                        && distance < best_distance
                    {
                        best = candidate;
                        best_distance = distance;
                    }
                }
            }
            best
        }
    }
}

pub fn multinomial_kernel(items: &[i64]) -> Result<f64, WorksheetErrorCode> {
    let positive_items = multinomial_positive_items(items);
    let exact_candidate = multinomial_exact_candidate(&positive_items);
    if !exact_candidate.is_finite() {
        return Ok(exact_candidate);
    }
    Ok(multinomial_publication_candidate(
        &positive_items,
        exact_candidate,
    ))
}