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