rand/distr/
uniform_float.rs

1// Copyright 2018-2020 Developers of the Rand project.
2// Copyright 2017 The Rust Project Developers.
3//
4// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
5// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
6// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
7// option. This file may not be copied, modified, or distributed
8// except according to those terms.
9
10//! `UniformFloat` implementation
11
12use super::{Error, SampleBorrow, SampleUniform, UniformSampler};
13use crate::distr::float::IntoFloat;
14use crate::distr::utils::{BoolAsSIMD, FloatAsSIMD, FloatSIMDUtils, IntAsSIMD};
15use crate::{Rng, RngExt};
16
17#[cfg(feature = "simd_support")]
18use core::simd::prelude::*;
19
20#[cfg(feature = "serde")]
21use serde::{Deserialize, Serialize};
22
23/// The back-end implementing [`UniformSampler`] for floating-point types.
24///
25/// Unless you are implementing [`UniformSampler`] for your own type, this type
26/// should not be used directly, use [`Uniform`] instead.
27///
28/// # Implementation notes
29///
30/// `UniformFloat` implementations convert RNG output to a float in the range
31/// `[1, 2)` via transmutation, map this to `[0, 1)`, then scale and translate
32/// to the desired range. Values produced this way have what equals 23 bits of
33/// random digits for an `f32` and 52 for an `f64`.
34///
35/// # Bias and range errors
36///
37/// Bias may be expected within the least-significant bit of the significand.
38/// It is not guaranteed that exclusive limits of a range are respected; i.e.
39/// when sampling the range `[a, b)` it is not guaranteed that `b` is never
40/// sampled.
41///
42/// [`new`]: UniformSampler::new
43/// [`new_inclusive`]: UniformSampler::new_inclusive
44/// [`StandardUniform`]: crate::distr::StandardUniform
45/// [`Uniform`]: super::Uniform
46#[derive(Clone, Copy, Debug, PartialEq)]
47#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
48pub struct UniformFloat<X> {
49    low: X,
50    scale: X,
51}
52
53macro_rules! uniform_float_impl {
54    ($($meta:meta)?, $ty:ty, $uty:ident, $f_scalar:ident, $u_scalar:ident, $bits_to_discard:expr) => {
55        $(#[cfg($meta)])?
56        impl UniformFloat<$ty> {
57            /// Construct, reducing `scale` as required to ensure that rounding
58            /// can never yield values greater than `high`.
59            ///
60            /// Requirements: `low` and `high` must be finite with `low <= high`.
61            /// `scale` must be non-negative, may be positive infinity but must not be NaN.
62            ///
63            /// Note: though it may be tempting to use a variant of this method
64            /// to ensure that samples from `[low, high)` are always strictly
65            /// less than `high`, this approach may be very slow where
66            /// `scale.abs()` is much smaller than `high.abs()`
67            /// (example: `low=0.99999999997819644, high=1.`).
68            fn new_bounded(low: $ty, high: $ty, mut scale: $ty) -> Self {
69                let max_rand = <$ty>::splat(1.0 as $f_scalar - $f_scalar::EPSILON);
70
71                loop {
72                    let mask = (scale * max_rand + low).le_mask(high);
73                    if mask.all() {
74                        break;
75                    }
76                    scale = scale.decrease_masked(!mask);
77                }
78
79                debug_assert!(<$ty>::splat(0.0).all_le(scale));
80
81                UniformFloat { low, scale }
82            }
83        }
84
85        $(#[cfg($meta)])?
86        impl SampleUniform for $ty {
87            type Sampler = UniformFloat<$ty>;
88        }
89
90        $(#[cfg($meta)])?
91        impl UniformSampler for UniformFloat<$ty> {
92            type X = $ty;
93
94            fn new<B1, B2>(low_b: B1, high_b: B2) -> Result<Self, Error>
95            where
96                B1: SampleBorrow<Self::X> + Sized,
97                B2: SampleBorrow<Self::X> + Sized,
98            {
99                let low = *low_b.borrow();
100                let high = *high_b.borrow();
101                #[cfg(debug_assertions)]
102                if !(low.all_finite()) || !(high.all_finite()) {
103                    return Err(Error::NonFinite);
104                }
105                if !(low.all_lt(high)) {
106                    return Err(Error::EmptyRange);
107                }
108
109                let scale = high - low;
110                if !(scale.all_finite()) {
111                    return Err(Error::NonFinite);
112                }
113
114                Ok(Self::new_bounded(low, high, scale))
115            }
116
117            fn new_inclusive<B1, B2>(low_b: B1, high_b: B2) -> Result<Self, Error>
118            where
119                B1: SampleBorrow<Self::X> + Sized,
120                B2: SampleBorrow<Self::X> + Sized,
121            {
122                let low = *low_b.borrow();
123                let high = *high_b.borrow();
124                #[cfg(debug_assertions)]
125                if !(low.all_finite()) || !(high.all_finite()) {
126                    return Err(Error::NonFinite);
127                }
128                if !low.all_le(high) {
129                    return Err(Error::EmptyRange);
130                }
131
132                let range = high - low;
133                if !range.all_finite() {
134                    return Err(Error::NonFinite);
135                }
136
137                let max_rand = <$ty>::splat(1.0 as $f_scalar - $f_scalar::EPSILON);
138                let scale = range / max_rand;
139
140                Ok(Self::new_bounded(low, high, scale))
141            }
142
143            fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::X {
144                // Generate a value in the range [1, 2)
145                let value1_2 = (rng.random::<$uty>() >> $uty::splat($bits_to_discard)).into_float_with_exponent(0);
146
147                // Get a value in the range [0, 1) to avoid overflow when multiplying by scale
148                let value0_1 = value1_2 - <$ty>::splat(1.0);
149
150                // We don't use `f64::mul_add`, because it is not available with
151                // `no_std`. Furthermore, it is slower for some targets (but
152                // faster for others). However, the order of multiplication and
153                // addition is important, because on some platforms (e.g. ARM)
154                // it will be optimized to a single (non-FMA) instruction.
155                value0_1 * self.scale + self.low
156            }
157
158            #[inline]
159            fn sample_single<R: Rng + ?Sized, B1, B2>(low_b: B1, high_b: B2, rng: &mut R) -> Result<Self::X, Error>
160            where
161                B1: SampleBorrow<Self::X> + Sized,
162                B2: SampleBorrow<Self::X> + Sized,
163            {
164                Self::sample_single_inclusive(low_b, high_b, rng)
165            }
166
167            #[inline]
168            fn sample_single_inclusive<R: Rng + ?Sized, B1, B2>(low_b: B1, high_b: B2, rng: &mut R) -> Result<Self::X, Error>
169            where
170                B1: SampleBorrow<Self::X> + Sized,
171                B2: SampleBorrow<Self::X> + Sized,
172            {
173                let low = *low_b.borrow();
174                let high = *high_b.borrow();
175                #[cfg(debug_assertions)]
176                if !low.all_finite() || !high.all_finite() {
177                    return Err(Error::NonFinite);
178                }
179                if !low.all_le(high) {
180                    return Err(Error::EmptyRange);
181                }
182                let scale = high - low;
183                if !scale.all_finite() {
184                    return Err(Error::NonFinite);
185                }
186
187                // Generate a value in the range [1, 2)
188                let value1_2 =
189                    (rng.random::<$uty>() >> $uty::splat($bits_to_discard)).into_float_with_exponent(0);
190
191                // Get a value in the range [0, 1) to avoid overflow when multiplying by scale
192                let value0_1 = value1_2 - <$ty>::splat(1.0);
193
194                // Doing multiply before addition allows some architectures
195                // to use a single instruction.
196                Ok(value0_1 * scale + low)
197            }
198        }
199    };
200}
201
202uniform_float_impl! { , f32, u32, f32, u32, 32 - 23 }
203uniform_float_impl! { , f64, u64, f64, u64, 64 - 52 }
204
205#[cfg(feature = "simd_support")]
206uniform_float_impl! { feature = "simd_support", f32x2, u32x2, f32, u32, 32 - 23 }
207#[cfg(feature = "simd_support")]
208uniform_float_impl! { feature = "simd_support", f32x4, u32x4, f32, u32, 32 - 23 }
209#[cfg(feature = "simd_support")]
210uniform_float_impl! { feature = "simd_support", f32x8, u32x8, f32, u32, 32 - 23 }
211#[cfg(feature = "simd_support")]
212uniform_float_impl! { feature = "simd_support", f32x16, u32x16, f32, u32, 32 - 23 }
213
214#[cfg(feature = "simd_support")]
215uniform_float_impl! { feature = "simd_support", f64x2, u64x2, f64, u64, 64 - 52 }
216#[cfg(feature = "simd_support")]
217uniform_float_impl! { feature = "simd_support", f64x4, u64x4, f64, u64, 64 - 52 }
218#[cfg(feature = "simd_support")]
219uniform_float_impl! { feature = "simd_support", f64x8, u64x8, f64, u64, 64 - 52 }
220
221#[cfg(test)]
222mod tests {
223    use super::*;
224    use crate::distr::{Uniform, utils::FloatSIMDScalarUtils};
225    use crate::test::{const_rng, step_rng};
226
227    #[test]
228    #[cfg_attr(miri, ignore)] // Miri is too slow
229    fn test_floats() {
230        let mut rng = crate::test::rng(252);
231        let mut zero_rng = const_rng(0);
232        let mut max_rng = const_rng(0xffff_ffff_ffff_ffff);
233        macro_rules! t {
234            ($ty:ty, $f_scalar:ident, $bits_shifted:expr) => {{
235                let v: &[($f_scalar, $f_scalar)] = &[
236                    (0.0, 100.0),
237                    (-1e35, -1e25),
238                    (1e-35, 1e-25),
239                    (-1e35, 1e35),
240                    (<$f_scalar>::from_bits(0), <$f_scalar>::from_bits(3)),
241                    (-<$f_scalar>::from_bits(10), -<$f_scalar>::from_bits(1)),
242                    (-<$f_scalar>::from_bits(5), 0.0),
243                    (-<$f_scalar>::from_bits(7), -0.0),
244                    (0.1 * $f_scalar::MAX, $f_scalar::MAX),
245                    (-$f_scalar::MAX * 0.2, $f_scalar::MAX * 0.7),
246                    (0.0, $f_scalar::MAX),
247                    (-$f_scalar::MAX, 0.0),
248                ];
249                for &(low_scalar, high_scalar) in v.iter() {
250                    for lane in 0..<$ty>::LEN {
251                        let low = <$ty>::splat(0.0 as $f_scalar).replace(lane, low_scalar);
252                        let high = <$ty>::splat(1.0 as $f_scalar).replace(lane, high_scalar);
253                        let my_uniform = Uniform::new(low, high).unwrap();
254                        let my_incl_uniform = Uniform::new_inclusive(low, high).unwrap();
255                        for _ in 0..100 {
256                            let v = rng.sample(my_uniform).extract_lane(lane);
257                            assert!(low_scalar <= v && v <= high_scalar);
258                            let v = rng.sample(my_incl_uniform).extract_lane(lane);
259                            assert!(low_scalar <= v && v <= high_scalar);
260                            let v =
261                                <$ty as SampleUniform>::Sampler::sample_single(low, high, &mut rng)
262                                    .unwrap()
263                                    .extract_lane(lane);
264                            assert!(low_scalar <= v && v <= high_scalar);
265                            let v = <$ty as SampleUniform>::Sampler::sample_single_inclusive(
266                                low, high, &mut rng,
267                            )
268                            .unwrap()
269                            .extract_lane(lane);
270                            assert!(low_scalar <= v && v <= high_scalar);
271                        }
272
273                        assert_eq!(
274                            rng.sample(Uniform::new_inclusive(low, low).unwrap())
275                                .extract_lane(lane),
276                            low_scalar
277                        );
278
279                        assert_eq!(zero_rng.sample(my_uniform).extract_lane(lane), low_scalar);
280                        assert_eq!(
281                            zero_rng.sample(my_incl_uniform).extract_lane(lane),
282                            low_scalar
283                        );
284                        assert_eq!(
285                            <$ty as SampleUniform>::Sampler::sample_single(
286                                low,
287                                high,
288                                &mut zero_rng
289                            )
290                            .unwrap()
291                            .extract_lane(lane),
292                            low_scalar
293                        );
294                        assert_eq!(
295                            <$ty as SampleUniform>::Sampler::sample_single_inclusive(
296                                low,
297                                high,
298                                &mut zero_rng
299                            )
300                            .unwrap()
301                            .extract_lane(lane),
302                            low_scalar
303                        );
304
305                        assert!(max_rng.sample(my_uniform).extract_lane(lane) <= high_scalar);
306                        assert!(max_rng.sample(my_incl_uniform).extract_lane(lane) <= high_scalar);
307                        // sample_single cannot cope with max_rng:
308                        // assert!(<$ty as SampleUniform>::Sampler
309                        //     ::sample_single(low, high, &mut max_rng).unwrap()
310                        //     .extract(lane) <= high_scalar);
311                        assert!(
312                            <$ty as SampleUniform>::Sampler::sample_single_inclusive(
313                                low,
314                                high,
315                                &mut max_rng
316                            )
317                            .unwrap()
318                            .extract_lane(lane)
319                                <= high_scalar
320                        );
321
322                        // Don't run this test for really tiny differences between high and low
323                        // since for those rounding might result in selecting high for a very
324                        // long time.
325                        if (high_scalar - low_scalar) > 0.0001 {
326                            let mut lowering_max_rng =
327                                step_rng(0xffff_ffff_ffff_ffff, (-1i64 << $bits_shifted) as u64);
328                            assert!(
329                                <$ty as SampleUniform>::Sampler::sample_single(
330                                    low,
331                                    high,
332                                    &mut lowering_max_rng
333                                )
334                                .unwrap()
335                                .extract_lane(lane)
336                                    <= high_scalar
337                            );
338                        }
339                    }
340                }
341
342                assert_eq!(
343                    rng.sample(Uniform::new_inclusive($f_scalar::MAX, $f_scalar::MAX).unwrap()),
344                    $f_scalar::MAX
345                );
346                assert_eq!(
347                    rng.sample(Uniform::new_inclusive(-$f_scalar::MAX, -$f_scalar::MAX).unwrap()),
348                    -$f_scalar::MAX
349                );
350            }};
351        }
352
353        t!(f32, f32, 32 - 23);
354        t!(f64, f64, 64 - 52);
355        #[cfg(feature = "simd_support")]
356        {
357            t!(f32x2, f32, 32 - 23);
358            t!(f32x4, f32, 32 - 23);
359            t!(f32x8, f32, 32 - 23);
360            t!(f32x16, f32, 32 - 23);
361            t!(f64x2, f64, 64 - 52);
362            t!(f64x4, f64, 64 - 52);
363            t!(f64x8, f64, 64 - 52);
364        }
365    }
366
367    #[test]
368    fn test_float_overflow() {
369        assert_eq!(Uniform::try_from(f64::MIN..f64::MAX), Err(Error::NonFinite));
370        assert_eq!(
371            Uniform::try_from(f64::MIN..=f64::MAX),
372            Err(Error::NonFinite)
373        );
374    }
375
376    #[test]
377    #[should_panic]
378    fn test_float_overflow_single() {
379        let mut rng = crate::test::rng(252);
380        rng.random_range(f64::MIN..f64::MAX);
381    }
382
383    #[test]
384    #[cfg(all(feature = "std", panic = "unwind"))]
385    fn test_float_assertions() {
386        use super::SampleUniform;
387        fn range<T: SampleUniform>(low: T, high: T) -> Result<T, Error> {
388            let mut rng = crate::test::rng(253);
389            T::Sampler::sample_single(low, high, &mut rng)
390        }
391
392        macro_rules! t {
393            ($ty:ident, $f_scalar:ident) => {{
394                let v: &[($f_scalar, $f_scalar)] = &[
395                    ($f_scalar::NAN, 0.0),
396                    (1.0, $f_scalar::NAN),
397                    ($f_scalar::NAN, $f_scalar::NAN),
398                    (1.0, 0.5),
399                    ($f_scalar::MAX, -$f_scalar::MAX),
400                    ($f_scalar::INFINITY, $f_scalar::INFINITY),
401                    ($f_scalar::NEG_INFINITY, $f_scalar::NEG_INFINITY),
402                    ($f_scalar::NEG_INFINITY, 5.0),
403                    (5.0, $f_scalar::INFINITY),
404                    ($f_scalar::NAN, $f_scalar::INFINITY),
405                    ($f_scalar::NEG_INFINITY, $f_scalar::NAN),
406                    ($f_scalar::NEG_INFINITY, $f_scalar::INFINITY),
407                ];
408                for &(low_scalar, high_scalar) in v.iter() {
409                    for lane in 0..<$ty>::LEN {
410                        let low = <$ty>::splat(0.0 as $f_scalar).replace(lane, low_scalar);
411                        let high = <$ty>::splat(1.0 as $f_scalar).replace(lane, high_scalar);
412                        assert!(range(low, high).is_err());
413                        assert!(Uniform::new(low, high).is_err());
414                        assert!(Uniform::new_inclusive(low, high).is_err());
415                        assert!(Uniform::new(low, low).is_err());
416                    }
417                }
418            }};
419        }
420
421        t!(f32, f32);
422        t!(f64, f64);
423        #[cfg(feature = "simd_support")]
424        {
425            t!(f32x2, f32);
426            t!(f32x4, f32);
427            t!(f32x8, f32);
428            t!(f32x16, f32);
429            t!(f64x2, f64);
430            t!(f64x4, f64);
431            t!(f64x8, f64);
432        }
433    }
434
435    #[test]
436    fn test_uniform_from_std_range() {
437        let r = Uniform::try_from(2.0f64..7.0).unwrap();
438        assert_eq!(r.0.low, 2.0);
439        assert_eq!(r.0.scale, 5.0);
440    }
441
442    #[test]
443    fn test_uniform_from_std_range_bad_limits() {
444        #![allow(clippy::reversed_empty_ranges)]
445        assert!(Uniform::try_from(100.0..10.0).is_err());
446        assert!(Uniform::try_from(100.0..100.0).is_err());
447    }
448
449    #[test]
450    fn test_uniform_from_std_range_inclusive() {
451        let r = Uniform::try_from(2.0f64..=7.0).unwrap();
452        assert_eq!(r.0.low, 2.0);
453        assert!(r.0.scale > 5.0);
454        assert!(r.0.scale < 5.0 + 1e-14);
455    }
456
457    #[test]
458    fn test_uniform_from_std_range_inclusive_bad_limits() {
459        #![allow(clippy::reversed_empty_ranges)]
460        assert!(Uniform::try_from(100.0..=10.0).is_err());
461        assert!(Uniform::try_from(100.0..=99.0).is_err());
462    }
463}
rand/distr/uniform_float.rs

rand/distr/
uniform_float.rs
