]>
Commit | Line | Data |
---|---|---|
0731742a | 1 | // Copyright 2018 Developers of the Rand project. |
b7449926 XL |
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 | ||
416331ca | 11 | use crate::distributions::utils::FloatSIMDUtils; |
dfeec247 XL |
12 | use crate::distributions::{Distribution, Standard}; |
13 | use crate::Rng; | |
14 | use core::mem; | |
15 | #[cfg(feature = "simd_support")] use packed_simd::*; | |
b7449926 XL |
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` | |
dfeec247 | 21 | /// the 24 most significant random bits of a `u32` are used and for `f64` the |
b7449926 XL |
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 | /// | |
416331ca XL |
38 | /// [`Standard`]: crate::distributions::Standard |
39 | /// [`Open01`]: crate::distributions::Open01 | |
40 | /// [`Uniform`]: crate::distributions::uniform::Uniform | |
b7449926 XL |
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` | |
dfeec247 | 48 | /// the 23 most significant random bits of an `u32` are used, for `f64` 52 from |
b7449926 XL |
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 | /// | |
416331ca XL |
64 | /// [`Standard`]: crate::distributions::Standard |
65 | /// [`OpenClosed01`]: crate::distributions::OpenClosed01 | |
66 | /// [`Uniform`]: crate::distributions::uniform::Uniform | |
b7449926 XL |
67 | #[derive(Clone, Copy, Debug)] |
68 | pub struct Open01; | |
69 | ||
70 | ||
416331ca XL |
71 | // This trait is needed by both this lib and rand_distr hence is a hidden export |
72 | #[doc(hidden)] | |
73 | pub trait IntoFloat { | |
b7449926 XL |
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). | |
84 | fn into_float_with_exponent(self, exponent: i32) -> Self::F; | |
85 | } | |
86 | ||
87 | macro_rules! float_impls { | |
0731742a XL |
88 | ($ty:ident, $uty:ident, $f_scalar:ident, $u_scalar:ty, |
89 | $fraction_bits:expr, $exponent_bias:expr) => { | |
b7449926 XL |
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 | |
0731742a XL |
95 | let exponent_bits: $u_scalar = |
96 | (($exponent_bias + exponent) as $u_scalar) << $fraction_bits; | |
dfeec247 | 97 | $ty::from_bits(self | exponent_bits) |
b7449926 XL |
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. | |
0731742a | 106 | let float_size = mem::size_of::<$f_scalar>() as u32 * 8; |
b7449926 | 107 | let precision = $fraction_bits + 1; |
0731742a | 108 | let scale = 1.0 / ((1 as $u_scalar << precision) as $f_scalar); |
b7449926 XL |
109 | |
110 | let value: $uty = rng.gen(); | |
0731742a XL |
111 | let value = value >> (float_size - precision); |
112 | scale * $ty::cast_from_int(value) | |
b7449926 XL |
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. | |
0731742a | 121 | let float_size = mem::size_of::<$f_scalar>() as u32 * 8; |
b7449926 | 122 | let precision = $fraction_bits + 1; |
0731742a | 123 | let scale = 1.0 / ((1 as $u_scalar << precision) as $f_scalar); |
b7449926 XL |
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: | |
0731742a | 128 | scale * $ty::cast_from_int(value + 1) |
b7449926 XL |
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. | |
0731742a XL |
137 | use core::$f_scalar::EPSILON; |
138 | let float_size = mem::size_of::<$f_scalar>() as u32 * 8; | |
b7449926 XL |
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 | } | |
0731742a XL |
147 | |
148 | float_impls! { f32, u32, f32, u32, 23, 127 } | |
149 | float_impls! { f64, u64, f64, u64, 52, 1023 } | |
150 | ||
dfeec247 | 151 | #[cfg(feature = "simd_support")] |
0731742a | 152 | float_impls! { f32x2, u32x2, f32, u32, 23, 127 } |
dfeec247 | 153 | #[cfg(feature = "simd_support")] |
0731742a | 154 | float_impls! { f32x4, u32x4, f32, u32, 23, 127 } |
dfeec247 | 155 | #[cfg(feature = "simd_support")] |
0731742a | 156 | float_impls! { f32x8, u32x8, f32, u32, 23, 127 } |
dfeec247 | 157 | #[cfg(feature = "simd_support")] |
0731742a XL |
158 | float_impls! { f32x16, u32x16, f32, u32, 23, 127 } |
159 | ||
dfeec247 | 160 | #[cfg(feature = "simd_support")] |
0731742a | 161 | float_impls! { f64x2, u64x2, f64, u64, 52, 1023 } |
dfeec247 | 162 | #[cfg(feature = "simd_support")] |
0731742a | 163 | float_impls! { f64x4, u64x4, f64, u64, 52, 1023 } |
dfeec247 | 164 | #[cfg(feature = "simd_support")] |
0731742a | 165 | float_impls! { f64x8, u64x8, f64, u64, 52, 1023 } |
b7449926 XL |
166 | |
167 | ||
168 | #[cfg(test)] | |
169 | mod tests { | |
dfeec247 | 170 | use super::*; |
416331ca | 171 | use crate::rngs::mock::StepRng; |
b7449926 XL |
172 | |
173 | const EPSILON32: f32 = ::core::f32::EPSILON; | |
174 | const EPSILON64: f64 = ::core::f64::EPSILON; | |
175 | ||
0731742a XL |
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); | |
b7449926 | 187 | |
0731742a XL |
188 | // OpenClosed01 |
189 | let mut zeros = StepRng::new(0, 0); | |
dfeec247 | 190 | assert_eq!(zeros.sample::<$ty, _>(OpenClosed01), 0.0 + $EPSILON / 2.0); |
0731742a XL |
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); | |
b7449926 | 195 | |
0731742a XL |
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 | } | |
dfeec247 | 204 | }; |
b7449926 | 205 | } |
0731742a | 206 | test_f32! { f32_edge_cases, f32, 0.0, EPSILON32 } |
dfeec247 | 207 | #[cfg(feature = "simd_support")] |
0731742a | 208 | test_f32! { f32x2_edge_cases, f32x2, f32x2::splat(0.0), f32x2::splat(EPSILON32) } |
dfeec247 | 209 | #[cfg(feature = "simd_support")] |
0731742a | 210 | test_f32! { f32x4_edge_cases, f32x4, f32x4::splat(0.0), f32x4::splat(EPSILON32) } |
dfeec247 | 211 | #[cfg(feature = "simd_support")] |
0731742a | 212 | test_f32! { f32x8_edge_cases, f32x8, f32x8::splat(0.0), f32x8::splat(EPSILON32) } |
dfeec247 | 213 | #[cfg(feature = "simd_support")] |
0731742a | 214 | test_f32! { f32x16_edge_cases, f32x16, f32x16::splat(0.0), f32x16::splat(EPSILON32) } |
b7449926 | 215 | |
0731742a XL |
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); | |
b7449926 | 227 | |
0731742a XL |
228 | // OpenClosed01 |
229 | let mut zeros = StepRng::new(0, 0); | |
dfeec247 | 230 | assert_eq!(zeros.sample::<$ty, _>(OpenClosed01), 0.0 + $EPSILON / 2.0); |
0731742a XL |
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); | |
b7449926 | 235 | |
0731742a XL |
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 | } | |
dfeec247 | 244 | }; |
b7449926 | 245 | } |
0731742a | 246 | test_f64! { f64_edge_cases, f64, 0.0, EPSILON64 } |
dfeec247 | 247 | #[cfg(feature = "simd_support")] |
0731742a | 248 | test_f64! { f64x2_edge_cases, f64x2, f64x2::splat(0.0), f64x2::splat(EPSILON64) } |
dfeec247 | 249 | #[cfg(feature = "simd_support")] |
0731742a | 250 | test_f64! { f64x4_edge_cases, f64x4, f64x4::splat(0.0), f64x4::splat(EPSILON64) } |
dfeec247 | 251 | #[cfg(feature = "simd_support")] |
0731742a | 252 | test_f64! { f64x8_edge_cases, f64x8, f64x8::splat(0.0), f64x8::splat(EPSILON64) } |
dfeec247 XL |
253 | |
254 | #[test] | |
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 | } | |
b7449926 | 307 | } |