use crate::ziggurat_tables;
use rand::distributions::hidden_export::IntoFloat;
use rand::Rng;
pub(crate) fn log_gamma<F: num_traits::Float>(x: F) -> F {
let coefficients: [F; 6] = [
F::from(76.18009172947146).unwrap(),
F::from(-86.50532032941677).unwrap(),
F::from(24.01409824083091).unwrap(),
F::from(-1.231739572450155).unwrap(),
F::from(0.1208650973866179e-2).unwrap(),
F::from(-0.5395239384953e-5).unwrap(),
];
let tmp = x + F::from(5.5).unwrap();
let log = (x + F::from(0.5).unwrap()) * tmp.ln() - tmp;
let mut a = F::from(1.000000000190015).unwrap();
let mut denom = x;
for &coeff in &coefficients {
denom = denom + F::one();
a = a + (coeff / denom);
}
log + (F::from(2.5066282746310005).unwrap() * a / x).ln()
}
#[inline(always)]
pub(crate) fn ziggurat<R: Rng + ?Sized, P, Z>(
rng: &mut R,
symmetric: bool,
x_tab: ziggurat_tables::ZigTable,
f_tab: ziggurat_tables::ZigTable,
mut pdf: P,
mut zero_case: Z
) -> f64
where
P: FnMut(f64) -> f64,
Z: FnMut(&mut R, f64) -> f64,
{
loop {
let bits = rng.next_u64();
let i = bits as usize & 0xff;
let u = if symmetric {
(bits >> 12).into_float_with_exponent(1) - 3.0
} else {
(bits >> 12).into_float_with_exponent(0) - (1.0 - core::f64::EPSILON / 2.0)
};
let x = u * x_tab[i];
let test_x = if symmetric { x.abs() } else { x };
if test_x < x_tab[i + 1] {
return x;
}
if i == 0 {
return zero_case(rng, u);
}
if f_tab[i + 1] + (f_tab[i] - f_tab[i + 1]) * rng.gen::<f64>() < pdf(x) {
return x;
}
}
}