|
| 1 | +use core::ops::{Add, Bound, Div, Mul, Neg, RangeBounds, Sub}; |
| 2 | + |
| 3 | +use crate::BaseRng; |
| 4 | + |
| 5 | +pub(super) fn f32(rng: &mut impl BaseRng, range: impl RangeBounds<f32>) -> f32 { |
| 6 | + float_range_impl(rng, range) |
| 7 | +} |
| 8 | + |
| 9 | +pub(super) fn f64(rng: &mut impl BaseRng, range: impl RangeBounds<f64>) -> f64 { |
| 10 | + float_range_impl(rng, range) |
| 11 | +} |
| 12 | + |
| 13 | +trait FloatExt: |
| 14 | + Add<Self, Output = Self> |
| 15 | + + Sub<Self, Output = Self> |
| 16 | + + Mul<Self, Output = Self> |
| 17 | + + Div<Self, Output = Self> |
| 18 | + + Neg<Output = Self> |
| 19 | + + Copy |
| 20 | + + Sized |
| 21 | +{ |
| 22 | + const MIN: Self; |
| 23 | + const MAX: Self; |
| 24 | + const HALF: Self; |
| 25 | + |
| 26 | + fn is_finite(self) -> bool; |
| 27 | + /// Generate a random float in [0, 1) range. |
| 28 | + fn gen_close_01_open(rng: &mut impl BaseRng) -> Self; |
| 29 | + /// Generate a random float in (0, 1] range. |
| 30 | + fn gen_open_01_close(rng: &mut impl BaseRng) -> Self; |
| 31 | + /// Generate a random float in (0, 1) range. |
| 32 | + fn gen_open_01_open(rng: &mut impl BaseRng) -> Self; |
| 33 | + /// Get the maximum float that can be generated in [0, 1) range (i.e., the |
| 34 | + /// float just below 1). |
| 35 | + fn max_rand() -> Self; |
| 36 | +} |
| 37 | + |
| 38 | +macro_rules! impl_float_ext { |
| 39 | + ($float:ident, $max_rand_bits:literal) => { |
| 40 | + impl FloatExt for $float { |
| 41 | + const MIN: Self = $float::MIN; |
| 42 | + const MAX: Self = $float::MAX; |
| 43 | + const HALF: Self = 0.5; |
| 44 | + |
| 45 | + #[inline] |
| 46 | + fn is_finite(self) -> bool { |
| 47 | + $float::is_finite(self) |
| 48 | + } |
| 49 | + #[inline] |
| 50 | + fn gen_close_01_open(rng: &mut impl BaseRng) -> Self { |
| 51 | + rng.$float() |
| 52 | + } |
| 53 | + #[inline] |
| 54 | + fn gen_open_01_close(rng: &mut impl BaseRng) -> Self { |
| 55 | + 1.0 - rng.$float() |
| 56 | + } |
| 57 | + #[inline] |
| 58 | + fn gen_open_01_open(rng: &mut impl BaseRng) -> Self { |
| 59 | + loop { |
| 60 | + let r = rng.$float(); |
| 61 | + |
| 62 | + if r != 0.0 { |
| 63 | + return r; |
| 64 | + } |
| 65 | + } |
| 66 | + } |
| 67 | + #[inline] |
| 68 | + fn max_rand() -> Self { |
| 69 | + <$float>::from_bits($max_rand_bits) |
| 70 | + } |
| 71 | + } |
| 72 | + }; |
| 73 | +} |
| 74 | + |
| 75 | +// Max rand constant is conceptually 0.9999... in the precision of the |
| 76 | +// particular float type. It can be determined in the following steps: |
| 77 | +// |
| 78 | +// 1. Fill the fractional part of the float with 1s (in binary). |
| 79 | +// 2. Move the number to have exponent 0 (i.e., the exponent is equal to the |
| 80 | +// exponent bias). |
| 81 | +// 3. Since the float representation uses implicit leading 1, subtract 1.0 |
| 82 | +// from the number constructed in the previous steps. |
| 83 | +// |
| 84 | +// This is a code snippet for f32: |
| 85 | +// |
| 86 | +// ``` |
| 87 | +// let fraction_bits = 23; // Significand bits without the implicit leading 1. |
| 88 | +// let exponent_bias = 127; |
| 89 | +// let discard_bits = u32::BITS - fraction_bits; |
| 90 | +// let max_rand = f32::from_bits(((u32::MAX >> discard_bits) | (exponent_bias << fraction_bits))) - 1.0; |
| 91 | +// ``` |
| 92 | + |
| 93 | +impl_float_ext!(f32, 0x3f7ffffe); |
| 94 | +impl_float_ext!(f64, 0x3feffffffffffffe); |
| 95 | + |
| 96 | +/// Indication whether a range is exclusive, inclusive on either side or |
| 97 | +/// inclusive on both sides. `Bound::Unbounded` is treated as inclusive, because |
| 98 | +/// we use MIN and MAX constants for such bounds. |
| 99 | +enum Inclusive { |
| 100 | + None, |
| 101 | + Left, |
| 102 | + Right, |
| 103 | + Both, |
| 104 | +} |
| 105 | + |
| 106 | +impl Inclusive { |
| 107 | + fn from_bounds<T>(range: impl RangeBounds<T>) -> Self { |
| 108 | + match (range.start_bound(), range.end_bound()) { |
| 109 | + (Bound::Excluded(_), Bound::Excluded(_)) => Self::None, |
| 110 | + (Bound::Included(_), Bound::Excluded(_)) | (Bound::Unbounded, Bound::Excluded(_)) => { |
| 111 | + Self::Left |
| 112 | + } |
| 113 | + (Bound::Excluded(_), Bound::Included(_)) | (Bound::Excluded(_), Bound::Unbounded) => { |
| 114 | + Self::Right |
| 115 | + } |
| 116 | + (Bound::Included(_), Bound::Included(_)) |
| 117 | + | (Bound::Included(_), Bound::Unbounded) |
| 118 | + | (Bound::Unbounded, Bound::Included(_)) |
| 119 | + | (Bound::Unbounded, Bound::Unbounded) => Self::Both, |
| 120 | + } |
| 121 | + } |
| 122 | +} |
| 123 | + |
| 124 | +fn float_range_impl<T: FloatExt>(rng: &mut impl BaseRng, range: impl RangeBounds<T>) -> T { |
| 125 | + let low = match range.start_bound() { |
| 126 | + Bound::Included(&low) | Bound::Excluded(&low) => low, |
| 127 | + Bound::Unbounded => T::MIN, |
| 128 | + }; |
| 129 | + |
| 130 | + let high = match range.end_bound() { |
| 131 | + Bound::Included(&high) | Bound::Excluded(&high) => high, |
| 132 | + Bound::Unbounded => T::MAX, |
| 133 | + }; |
| 134 | + |
| 135 | + let inclusive = Inclusive::from_bounds(range); |
| 136 | + |
| 137 | + // Our generator is able to generate floats with one or both sides of the |
| 138 | + // range open. However, it can't generate a float from the range closed on |
| 139 | + // both sides. For this case, we divide the scale by maximum random number |
| 140 | + // which "stretches" the range to include both sides. This is the approach |
| 141 | + // used in rand crate: |
| 142 | + // https://github.com/rust-random/rand/blob/f3dd0b885c4597b9617ca79987a0dd899ab29fcb/src/distributions/uniform.rs#L953 |
| 143 | + let scale = match inclusive { |
| 144 | + Inclusive::None | Inclusive::Left | Inclusive::Right => high - low, |
| 145 | + Inclusive::Both => (high - low) / T::max_rand(), |
| 146 | + }; |
| 147 | + |
| 148 | + if scale.is_finite() { |
| 149 | + // Generate a random number between 0 and 1, where the bounds are |
| 150 | + // included based on the desired range inclusiveness. |
| 151 | + let r = match inclusive { |
| 152 | + Inclusive::None => T::gen_open_01_open(rng), |
| 153 | + Inclusive::Right => T::gen_open_01_close(rng), |
| 154 | + Inclusive::Left => T::gen_close_01_open(rng), |
| 155 | + // Inclusiveness on both sides is achieved by stretching the scale |
| 156 | + // above. |
| 157 | + Inclusive::Both => T::gen_close_01_open(rng), |
| 158 | + }; |
| 159 | + |
| 160 | + r * scale + low |
| 161 | + } else { |
| 162 | + // Scale not being finite means that the range is wider than the float |
| 163 | + // type can represent (or that at least one side is not finite). In such |
| 164 | + // case, we need to fall back into the following technique which does a |
| 165 | + // bit more work but can handle such ranges. Source: |
| 166 | + // https://medium.com/analytics-vidhya/random-floats-in-any-range-9b40d30b637b |
| 167 | + let high_half = T::HALF * high; |
| 168 | + let low_half = T::HALF * low; |
| 169 | + let mid_point = high_half + low_half; |
| 170 | + |
| 171 | + // Decide if we generate the value to the right or left from the middle |
| 172 | + // point. We always want to have a chance that the middle point is |
| 173 | + // sampled, so we can't use the (0, 1] trick with one-side inclusive |
| 174 | + // ranges. That is why we stretch those in appropriate cases. |
| 175 | + let (r, stretch) = if rng.bool() { |
| 176 | + let stretch = match inclusive { |
| 177 | + Inclusive::None | Inclusive::Left => false, |
| 178 | + Inclusive::Right | Inclusive::Both => true, |
| 179 | + }; |
| 180 | + let r = T::gen_close_01_open(rng); |
| 181 | + (r, stretch) |
| 182 | + } else { |
| 183 | + let stretch = match inclusive { |
| 184 | + Inclusive::None | Inclusive::Right => false, |
| 185 | + Inclusive::Left | Inclusive::Both => true, |
| 186 | + }; |
| 187 | + let r = -T::gen_close_01_open(rng); |
| 188 | + (r, stretch) |
| 189 | + }; |
| 190 | + |
| 191 | + let half_scale = if stretch { |
| 192 | + let half_scale = (high_half - low_half) / T::max_rand(); |
| 193 | + if half_scale.is_finite() { |
| 194 | + half_scale |
| 195 | + } else { |
| 196 | + // If the range is so extreme that it can't be stretched, |
| 197 | + // use the standard scale. |
| 198 | + high_half - low_half |
| 199 | + } |
| 200 | + } else { |
| 201 | + high_half - low_half |
| 202 | + }; |
| 203 | + |
| 204 | + r * half_scale + mid_point |
| 205 | + } |
| 206 | +} |
| 207 | + |
| 208 | +#[cfg(test)] |
| 209 | +mod tests { |
| 210 | + use fastrand::Rng; |
| 211 | + |
| 212 | + use super::*; |
| 213 | + |
| 214 | + #[test] |
| 215 | + fn f32_range_in_bounds() { |
| 216 | + let mut rng = Rng::new(); |
| 217 | + |
| 218 | + let range = -2.0..2.0; |
| 219 | + for _ in 0..10000 { |
| 220 | + assert!(range.contains(&float_range_impl(&mut rng, range.clone()))); |
| 221 | + } |
| 222 | + } |
| 223 | + |
| 224 | + #[test] |
| 225 | + fn f32_range_wide_range_in_bounds() { |
| 226 | + let mut rng = Rng::new(); |
| 227 | + |
| 228 | + let range = f32::MIN..f32::MAX; |
| 229 | + for _ in 0..10000 { |
| 230 | + assert!(range.contains(&float_range_impl(&mut rng, range.clone()))); |
| 231 | + } |
| 232 | + } |
| 233 | + |
| 234 | + #[test] |
| 235 | + fn f32_range_unbounded_finite() { |
| 236 | + let mut rng = Rng::new(); |
| 237 | + |
| 238 | + let range = ..; |
| 239 | + for _ in 0..10000 { |
| 240 | + assert!(&float_range_impl::<f32>(&mut rng, range).is_finite()); |
| 241 | + } |
| 242 | + } |
| 243 | +} |
0 commit comments