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 //! Basic floating-point number distributions
10 
11 use crate::distributions::utils::FloatSIMDUtils;
12 use crate::distributions::{Distribution, Standard};
13 use crate::Rng;
14 use core::mem;
15 #[cfg(feature = "simd_support")] use packed_simd::*;
16 
17 /// A distribution to sample floating point numbers uniformly in the half-open
18 /// interval `(0, 1]`, i.e. including 1 but not 0.
19 ///
20 /// All values that can be generated are of the form `n * ε/2`. For `f32`
21 /// the 24 most significant random bits of a `u32` are used and for `f64` the
22 /// 53 most significant bits of a `u64` are used. The conversion uses the
23 /// multiplicative method.
24 ///
25 /// See also: [`Standard`] which samples from `[0, 1)`, [`Open01`]
26 /// which samples from `(0, 1)` and [`Uniform`] which samples from arbitrary
27 /// ranges.
28 ///
29 /// # Example
30 /// ```
31 /// use rand::{thread_rng, Rng};
32 /// use rand::distributions::OpenClosed01;
33 ///
34 /// let val: f32 = thread_rng().sample(OpenClosed01);
35 /// println!("f32 from (0, 1): {}", val);
36 /// ```
37 ///
38 /// [`Standard`]: crate::distributions::Standard
39 /// [`Open01`]: crate::distributions::Open01
40 /// [`Uniform`]: crate::distributions::uniform::Uniform
41 #[derive(Clone, Copy, Debug)]
42 pub struct OpenClosed01;
43 
44 /// A distribution to sample floating point numbers uniformly in the open
45 /// interval `(0, 1)`, i.e. not including either endpoint.
46 ///
47 /// All values that can be generated are of the form `n * ε + ε/2`. For `f32`
48 /// the 23 most significant random bits of an `u32` are used, for `f64` 52 from
49 /// an `u64`. The conversion uses a transmute-based method.
50 ///
51 /// See also: [`Standard`] which samples from `[0, 1)`, [`OpenClosed01`]
52 /// which samples from `(0, 1]` and [`Uniform`] which samples from arbitrary
53 /// ranges.
54 ///
55 /// # Example
56 /// ```
57 /// use rand::{thread_rng, Rng};
58 /// use rand::distributions::Open01;
59 ///
60 /// let val: f32 = thread_rng().sample(Open01);
61 /// println!("f32 from (0, 1): {}", val);
62 /// ```
63 ///
64 /// [`Standard`]: crate::distributions::Standard
65 /// [`OpenClosed01`]: crate::distributions::OpenClosed01
66 /// [`Uniform`]: crate::distributions::uniform::Uniform
67 #[derive(Clone, Copy, Debug)]
68 pub struct Open01;
69 
70 
71 // This trait is needed by both this lib and rand_distr hence is a hidden export
72 #[doc(hidden)]
73 pub trait IntoFloat {
74     type F;
75 
76     /// Helper method to combine the fraction and a contant exponent into a
77     /// float.
78     ///
79     /// Only the least significant bits of `self` may be set, 23 for `f32` and
80     /// 52 for `f64`.
81     /// The resulting value will fall in a range that depends on the exponent.
82     /// As an example the range with exponent 0 will be
83     /// [2<sup>0</sup>..2<sup>1</sup>), which is [1..2).
into_float_with_exponent(self, exponent: i32) -> Self::F84     fn into_float_with_exponent(self, exponent: i32) -> Self::F;
85 }
86 
87 macro_rules! float_impls {
88     ($ty:ident, $uty:ident, $f_scalar:ident, $u_scalar:ty,
89      $fraction_bits:expr, $exponent_bias:expr) => {
90         impl IntoFloat for $uty {
91             type F = $ty;
92             #[inline(always)]
93             fn into_float_with_exponent(self, exponent: i32) -> $ty {
94                 // The exponent is encoded using an offset-binary representation
95                 let exponent_bits: $u_scalar =
96                     (($exponent_bias + exponent) as $u_scalar) << $fraction_bits;
97                 $ty::from_bits(self | exponent_bits)
98             }
99         }
100 
101         impl Distribution<$ty> for Standard {
102             fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> $ty {
103                 // Multiply-based method; 24/53 random bits; [0, 1) interval.
104                 // We use the most significant bits because for simple RNGs
105                 // those are usually more random.
106                 let float_size = mem::size_of::<$f_scalar>() as u32 * 8;
107                 let precision = $fraction_bits + 1;
108                 let scale = 1.0 / ((1 as $u_scalar << precision) as $f_scalar);
109 
110                 let value: $uty = rng.gen();
111                 let value = value >> (float_size - precision);
112                 scale * $ty::cast_from_int(value)
113             }
114         }
115 
116         impl Distribution<$ty> for OpenClosed01 {
117             fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> $ty {
118                 // Multiply-based method; 24/53 random bits; (0, 1] interval.
119                 // We use the most significant bits because for simple RNGs
120                 // those are usually more random.
121                 let float_size = mem::size_of::<$f_scalar>() as u32 * 8;
122                 let precision = $fraction_bits + 1;
123                 let scale = 1.0 / ((1 as $u_scalar << precision) as $f_scalar);
124 
125                 let value: $uty = rng.gen();
126                 let value = value >> (float_size - precision);
127                 // Add 1 to shift up; will not overflow because of right-shift:
128                 scale * $ty::cast_from_int(value + 1)
129             }
130         }
131 
132         impl Distribution<$ty> for Open01 {
133             fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> $ty {
134                 // Transmute-based method; 23/52 random bits; (0, 1) interval.
135                 // We use the most significant bits because for simple RNGs
136                 // those are usually more random.
137                 use core::$f_scalar::EPSILON;
138                 let float_size = mem::size_of::<$f_scalar>() as u32 * 8;
139 
140                 let value: $uty = rng.gen();
141                 let fraction = value >> (float_size - $fraction_bits);
142                 fraction.into_float_with_exponent(0) - (1.0 - EPSILON / 2.0)
143             }
144         }
145     }
146 }
147 
148 float_impls! { f32, u32, f32, u32, 23, 127 }
149 float_impls! { f64, u64, f64, u64, 52, 1023 }
150 
151 #[cfg(feature = "simd_support")]
152 float_impls! { f32x2, u32x2, f32, u32, 23, 127 }
153 #[cfg(feature = "simd_support")]
154 float_impls! { f32x4, u32x4, f32, u32, 23, 127 }
155 #[cfg(feature = "simd_support")]
156 float_impls! { f32x8, u32x8, f32, u32, 23, 127 }
157 #[cfg(feature = "simd_support")]
158 float_impls! { f32x16, u32x16, f32, u32, 23, 127 }
159 
160 #[cfg(feature = "simd_support")]
161 float_impls! { f64x2, u64x2, f64, u64, 52, 1023 }
162 #[cfg(feature = "simd_support")]
163 float_impls! { f64x4, u64x4, f64, u64, 52, 1023 }
164 #[cfg(feature = "simd_support")]
165 float_impls! { f64x8, u64x8, f64, u64, 52, 1023 }
166 
167 
168 #[cfg(test)]
169 mod tests {
170     use super::*;
171     use crate::rngs::mock::StepRng;
172 
173     const EPSILON32: f32 = ::core::f32::EPSILON;
174     const EPSILON64: f64 = ::core::f64::EPSILON;
175 
176     macro_rules! test_f32 {
177         ($fnn:ident, $ty:ident, $ZERO:expr, $EPSILON:expr) => {
178             #[test]
179             fn $fnn() {
180                 // Standard
181                 let mut zeros = StepRng::new(0, 0);
182                 assert_eq!(zeros.gen::<$ty>(), $ZERO);
183                 let mut one = StepRng::new(1 << 8 | 1 << (8 + 32), 0);
184                 assert_eq!(one.gen::<$ty>(), $EPSILON / 2.0);
185                 let mut max = StepRng::new(!0, 0);
186                 assert_eq!(max.gen::<$ty>(), 1.0 - $EPSILON / 2.0);
187 
188                 // OpenClosed01
189                 let mut zeros = StepRng::new(0, 0);
190                 assert_eq!(zeros.sample::<$ty, _>(OpenClosed01), 0.0 + $EPSILON / 2.0);
191                 let mut one = StepRng::new(1 << 8 | 1 << (8 + 32), 0);
192                 assert_eq!(one.sample::<$ty, _>(OpenClosed01), $EPSILON);
193                 let mut max = StepRng::new(!0, 0);
194                 assert_eq!(max.sample::<$ty, _>(OpenClosed01), $ZERO + 1.0);
195 
196                 // Open01
197                 let mut zeros = StepRng::new(0, 0);
198                 assert_eq!(zeros.sample::<$ty, _>(Open01), 0.0 + $EPSILON / 2.0);
199                 let mut one = StepRng::new(1 << 9 | 1 << (9 + 32), 0);
200                 assert_eq!(one.sample::<$ty, _>(Open01), $EPSILON / 2.0 * 3.0);
201                 let mut max = StepRng::new(!0, 0);
202                 assert_eq!(max.sample::<$ty, _>(Open01), 1.0 - $EPSILON / 2.0);
203             }
204         };
205     }
206     test_f32! { f32_edge_cases, f32, 0.0, EPSILON32 }
207     #[cfg(feature = "simd_support")]
208     test_f32! { f32x2_edge_cases, f32x2, f32x2::splat(0.0), f32x2::splat(EPSILON32) }
209     #[cfg(feature = "simd_support")]
210     test_f32! { f32x4_edge_cases, f32x4, f32x4::splat(0.0), f32x4::splat(EPSILON32) }
211     #[cfg(feature = "simd_support")]
212     test_f32! { f32x8_edge_cases, f32x8, f32x8::splat(0.0), f32x8::splat(EPSILON32) }
213     #[cfg(feature = "simd_support")]
214     test_f32! { f32x16_edge_cases, f32x16, f32x16::splat(0.0), f32x16::splat(EPSILON32) }
215 
216     macro_rules! test_f64 {
217         ($fnn:ident, $ty:ident, $ZERO:expr, $EPSILON:expr) => {
218             #[test]
219             fn $fnn() {
220                 // Standard
221                 let mut zeros = StepRng::new(0, 0);
222                 assert_eq!(zeros.gen::<$ty>(), $ZERO);
223                 let mut one = StepRng::new(1 << 11, 0);
224                 assert_eq!(one.gen::<$ty>(), $EPSILON / 2.0);
225                 let mut max = StepRng::new(!0, 0);
226                 assert_eq!(max.gen::<$ty>(), 1.0 - $EPSILON / 2.0);
227 
228                 // OpenClosed01
229                 let mut zeros = StepRng::new(0, 0);
230                 assert_eq!(zeros.sample::<$ty, _>(OpenClosed01), 0.0 + $EPSILON / 2.0);
231                 let mut one = StepRng::new(1 << 11, 0);
232                 assert_eq!(one.sample::<$ty, _>(OpenClosed01), $EPSILON);
233                 let mut max = StepRng::new(!0, 0);
234                 assert_eq!(max.sample::<$ty, _>(OpenClosed01), $ZERO + 1.0);
235 
236                 // Open01
237                 let mut zeros = StepRng::new(0, 0);
238                 assert_eq!(zeros.sample::<$ty, _>(Open01), 0.0 + $EPSILON / 2.0);
239                 let mut one = StepRng::new(1 << 12, 0);
240                 assert_eq!(one.sample::<$ty, _>(Open01), $EPSILON / 2.0 * 3.0);
241                 let mut max = StepRng::new(!0, 0);
242                 assert_eq!(max.sample::<$ty, _>(Open01), 1.0 - $EPSILON / 2.0);
243             }
244         };
245     }
246     test_f64! { f64_edge_cases, f64, 0.0, EPSILON64 }
247     #[cfg(feature = "simd_support")]
248     test_f64! { f64x2_edge_cases, f64x2, f64x2::splat(0.0), f64x2::splat(EPSILON64) }
249     #[cfg(feature = "simd_support")]
250     test_f64! { f64x4_edge_cases, f64x4, f64x4::splat(0.0), f64x4::splat(EPSILON64) }
251     #[cfg(feature = "simd_support")]
252     test_f64! { f64x8_edge_cases, f64x8, f64x8::splat(0.0), f64x8::splat(EPSILON64) }
253 
254     #[test]
value_stability()255     fn value_stability() {
256         fn test_samples<T: Copy + core::fmt::Debug + PartialEq, D: Distribution<T>>(
257             distr: &D, zero: T, expected: &[T],
258         ) {
259             let mut rng = crate::test::rng(0x6f44f5646c2a7334);
260             let mut buf = [zero; 3];
261             for x in &mut buf {
262                 *x = rng.sample(&distr);
263             }
264             assert_eq!(&buf, expected);
265         }
266 
267         test_samples(&Standard, 0f32, &[0.0035963655, 0.7346052, 0.09778172]);
268         test_samples(&Standard, 0f64, &[
269             0.7346051961657583,
270             0.20298547462974248,
271             0.8166436635290655,
272         ]);
273 
274         test_samples(&OpenClosed01, 0f32, &[0.003596425, 0.73460525, 0.09778178]);
275         test_samples(&OpenClosed01, 0f64, &[
276             0.7346051961657584,
277             0.2029854746297426,
278             0.8166436635290656,
279         ]);
280 
281         test_samples(&Open01, 0f32, &[0.0035963655, 0.73460525, 0.09778172]);
282         test_samples(&Open01, 0f64, &[
283             0.7346051961657584,
284             0.20298547462974248,
285             0.8166436635290656,
286         ]);
287 
288         #[cfg(feature = "simd_support")]
289         {
290             // We only test a sub-set of types here. Values are identical to
291             // non-SIMD types; we assume this pattern continues across all
292             // SIMD types.
293 
294             test_samples(&Standard, f32x2::new(0.0, 0.0), &[
295                 f32x2::new(0.0035963655, 0.7346052),
296                 f32x2::new(0.09778172, 0.20298547),
297                 f32x2::new(0.34296435, 0.81664366),
298             ]);
299 
300             test_samples(&Standard, f64x2::new(0.0, 0.0), &[
301                 f64x2::new(0.7346051961657583, 0.20298547462974248),
302                 f64x2::new(0.8166436635290655, 0.7423708925400552),
303                 f64x2::new(0.16387782224016323, 0.9087068770169618),
304             ]);
305         }
306     }
307 }
308