rand_distr/
utils.rs

1// Copyright 2018 Developers of the Rand project.
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
5// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
6// option. This file may not be copied, modified, or distributed
7// except according to those terms.
8
9//! Math helper functions
10
11use crate::ziggurat_tables;
12#[allow(unused_imports)]
13use num_traits::Float; // Used for `no_std` to get `f64::abs()` working before `rustc 1.84`
14use rand::distr::hidden_export::IntoFloat;
15use rand::Rng;
16
17/// Sample a random number using the Ziggurat method (specifically the
18/// ZIGNOR variant from Doornik 2005). Most of the arguments are
19/// directly from the paper:
20///
21/// * `rng`: source of randomness
22/// * `symmetric`: whether this is a symmetric distribution, or one-sided with P(x < 0) = 0.
23/// * `X`: the $x_i$ abscissae.
24/// * `F`: precomputed values of the PDF at the $x_i$, (i.e. $f(x_i)$)
25/// * `F_DIFF`: precomputed values of $f(x_i) - f(x_{i+1})$
26/// * `pdf`: the probability density function
27/// * `zero_case`: manual sampling from the tail when we chose the
28///    bottom box (i.e. i == 0)
29#[inline(always)] // Forced inlining improves the perf by 25-50%
30pub(crate) fn ziggurat<R: Rng + ?Sized, P, Z>(
31    rng: &mut R,
32    symmetric: bool,
33    x_tab: ziggurat_tables::ZigTable,
34    f_tab: ziggurat_tables::ZigTable,
35    mut pdf: P,
36    mut zero_case: Z,
37) -> f64
38where
39    P: FnMut(f64) -> f64,
40    Z: FnMut(&mut R, f64) -> f64,
41{
42    loop {
43        // As an optimisation we re-implement the conversion to a f64.
44        // From the remaining 12 most significant bits we use 8 to construct `i`.
45        // This saves us generating a whole extra random number, while the added
46        // precision of using 64 bits for f64 does not buy us much.
47        let bits = rng.next_u64();
48        let i = bits as usize & 0xff;
49
50        let u = if symmetric {
51            // Convert to a value in the range [2,4) and subtract to get [-1,1)
52            // We can't convert to an open range directly, that would require
53            // subtracting `3.0 - EPSILON`, which is not representable.
54            // It is possible with an extra step, but an open range does not
55            // seem necessary for the ziggurat algorithm anyway.
56            (bits >> 12).into_float_with_exponent(1) - 3.0
57        } else {
58            // Convert to a value in the range [1,2) and subtract to get (0,1)
59            (bits >> 12).into_float_with_exponent(0) - (1.0 - f64::EPSILON / 2.0)
60        };
61        let x = u * x_tab[i];
62
63        let test_x = if symmetric { x.abs() } else { x };
64
65        // algebraically equivalent to |u| < x_tab[i+1]/x_tab[i] (or u < x_tab[i+1]/x_tab[i])
66        if test_x < x_tab[i + 1] {
67            return x;
68        }
69        if i == 0 {
70            return zero_case(rng, u);
71        }
72        // algebraically equivalent to f1 + DRanU()*(f0 - f1) < 1
73        if f_tab[i + 1] + (f_tab[i] - f_tab[i + 1]) * rng.random::<f64>() < pdf(x) {
74            return x;
75        }
76    }
77}
rand_distr/utils.rs

rand_distr/
utils.rs
