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