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