// Copyright 2018 Developers of the Rand project.

//

// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or

// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license

// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your

// option. This file may not be copied, modified, or distributed

// except according to those terms.

//! Math helper functions

use crate::ziggurat_tables;

#[allow(unused_imports)]

use num_traits::Float; // Used for `no_std` to get `f64::abs()` working before `rustc 1.84`

use rand::distr::hidden_export::IntoFloat;

use rand::Rng;

/// Sample a random number using the Ziggurat method (specifically the

/// ZIGNOR variant from Doornik 2005). Most of the arguments are

/// directly from the paper:

///

/// * `rng`: source of randomness

/// * `symmetric`: whether this is a symmetric distribution, or one-sided with P(x < 0) = 0.

/// * `X`: the $x_i$ abscissae.

/// * `F`: precomputed values of the PDF at the $x_i$, (i.e. $f(x_i)$)

/// * `F_DIFF`: precomputed values of $f(x_i) - f(x_{i+1})$

/// * `pdf`: the probability density function

/// * `zero_case`: manual sampling from the tail when we chose the

///    bottom box (i.e. i == 0)

#[inline(always)] // Forced inlining improves the perf by 25-50%

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 {

        // As an optimisation we re-implement the conversion to a f64.

        // From the remaining 12 most significant bits we use 8 to construct `i`.

        // This saves us generating a whole extra random number, while the added

        // precision of using 64 bits for f64 does not buy us much.

        let bits = rng.next_u64();

        let i = bits as usize & 0xff;

        let u = if symmetric {

            // Convert to a value in the range [2,4) and subtract to get [-1,1)

            // We can't convert to an open range directly, that would require

            // subtracting `3.0 - EPSILON`, which is not representable.

            // It is possible with an extra step, but an open range does not

            // seem necessary for the ziggurat algorithm anyway.

            (bits >> 12).into_float_with_exponent(1) - 3.0

        } else {

            // Convert to a value in the range [1,2) and subtract to get (0,1)

            (bits >> 12).into_float_with_exponent(0) - (1.0 - f64::EPSILON / 2.0)

        };

        let x = u * x_tab[i];

        let test_x = if symmetric { x.abs() } else { x };

        // algebraically equivalent to |u| < x_tab[i+1]/x_tab[i] (or u < x_tab[i+1]/x_tab[i])

        if test_x < x_tab[i + 1] {

            return x;

        }

        if i == 0 {

            return zero_case(rng, u);

        }

        // algebraically equivalent to f1 + DRanU()*(f0 - f1) < 1

        if f_tab[i + 1] + (f_tab[i] - f_tab[i + 1]) * rng.random::<f64>() < pdf(x) {

            return x;

        }

    }

}

Unexecuted instantiation: rand_distr::utils::ziggurat::<rand::rngs::std::StdRng, <rand_distr::normal::StandardNormal as rand::distr::distribution::Distribution<f64>>::sample::pdf, <rand_distr::normal::StandardNormal as rand::distr::distribution::Distribution<f64>>::sample::zero_case<rand::rngs::std::StdRng>>

Unexecuted instantiation: rand_distr::utils::ziggurat::<_, _, _>