@@ -14,203 +14,86 @@ use core::mem;
14
14
use Rng ;
15
15
use distributions:: { Distribution , Uniform } ;
16
16
17
-
18
- /// A distribution to sample floating point numbers uniformly in the open
19
- /// interval `(0, 1)` (not including either endpoint).
20
- ///
21
- /// See also: [`Closed01`] for the closed `[0, 1]`; [`Uniform`] for the
22
- /// half-open `[0, 1)`.
23
- ///
24
- /// # Example
25
- /// ```rust
26
- /// use rand::{weak_rng, Rng};
27
- /// use rand::distributions::Open01;
28
- ///
29
- /// let val: f32 = weak_rng().sample(Open01);
30
- /// println!("f32 from (0,1): {}", val);
31
- /// ```
32
- ///
33
- /// [`Uniform`]: struct.Uniform.html
34
- /// [`Closed01`]: struct.Closed01.html
35
- #[ derive( Clone , Copy , Debug ) ]
36
- pub struct Open01 ;
37
-
38
- /// A distribution to sample floating point numbers uniformly in the closed
39
- /// interval `[0, 1]` (including both endpoints).
40
- ///
41
- /// See also: [`Open01`] for the open `(0, 1)`; [`Uniform`] for the half-open
42
- /// `[0, 1)`.
43
- ///
44
- /// # Example
45
- /// ```rust
46
- /// use rand::{weak_rng, Rng};
47
- /// use rand::distributions::Closed01;
48
- ///
49
- /// let val: f32 = weak_rng().sample(Closed01);
50
- /// println!("f32 from [0,1]: {}", val);
51
- /// ```
52
- ///
53
- /// [`Uniform`]: struct.Uniform.html
54
- /// [`Open01`]: struct.Open01.html
55
- #[ derive( Clone , Copy , Debug ) ]
56
- pub struct Closed01 ;
57
-
58
-
59
- // Return the next random f32 selected from the half-open
60
- // interval `[0, 1)`.
61
- //
62
- // This uses a technique described by Saito and Matsumoto at
63
- // MCQMC'08. Given that the IEEE floating point numbers are
64
- // uniformly distributed over [1,2), we generate a number in
65
- // this range and then offset it onto the range [0,1). Our
66
- // choice of bits (masking v. shifting) is arbitrary and
67
- // should be immaterial for high quality generators. For low
68
- // quality generators (ex. LCG), prefer bitshifting due to
69
- // correlation between sequential low order bits.
70
- //
71
- // See:
72
- // A PRNG specialized in double precision floating point numbers using
73
- // an affine transition
74
- //
75
- // * <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/dSFMT.pdf>
76
- // * <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/dSFMT-slide-e.pdf>
77
- impl Distribution < f32 > for Uniform {
78
- fn sample < R : Rng + ?Sized > ( & self , rng : & mut R ) -> f32 {
79
- const UPPER_MASK : u32 = 0x3F800000 ;
80
- const LOWER_MASK : u32 = 0x7FFFFF ;
81
- let tmp = UPPER_MASK | ( rng. next_u32 ( ) & LOWER_MASK ) ;
82
- let result: f32 = unsafe { mem:: transmute ( tmp) } ;
83
- result - 1.0
84
- }
85
- }
86
- impl Distribution < f64 > for Uniform {
87
- fn sample < R : Rng + ?Sized > ( & self , rng : & mut R ) -> f64 {
88
- const UPPER_MASK : u64 = 0x3FF0000000000000 ;
89
- const LOWER_MASK : u64 = 0xFFFFFFFFFFFFF ;
90
- let tmp = UPPER_MASK | ( rng. next_u64 ( ) & LOWER_MASK ) ;
91
- let result: f64 = unsafe { mem:: transmute ( tmp) } ;
92
- result - 1.0
93
- }
17
+ pub ( crate ) trait IntoFloat {
18
+ type F ;
19
+
20
+ /// Helper method to combine the fraction and a contant exponent into a
21
+ /// float.
22
+ ///
23
+ /// Only the least significant bits of `self` may be set, 23 for `f32` and
24
+ /// 52 for `f64`.
25
+ /// The resulting value will fall in a range that depends on the exponent.
26
+ /// As an example the range with exponent 0 will be
27
+ /// [2<sup>0</sup>..2<sup>1</sup>), which is [1..2).
28
+ #[ inline( always) ]
29
+ fn into_float_with_exponent ( self , exponent : i32 ) -> Self :: F ;
94
30
}
95
31
96
32
macro_rules! float_impls {
97
- ( $mod_name: ident, $ty: ty, $mantissa_bits: expr) => {
98
- mod $mod_name {
99
- use Rng ;
100
- use distributions:: { Distribution } ;
101
- use super :: { Open01 , Closed01 } ;
102
-
103
- const SCALE : $ty = ( 1u64 << $mantissa_bits) as $ty;
104
-
105
- impl Distribution <$ty> for Open01 {
106
- #[ inline]
107
- fn sample<R : Rng + ?Sized >( & self , rng: & mut R ) -> $ty {
108
- // add 0.5 * epsilon, so that smallest number is
109
- // greater than 0, and largest number is still
110
- // less than 1, specifically 1 - 0.5 * epsilon.
111
- let x: $ty = rng. gen ( ) ;
112
- x + 0.5 / SCALE
113
- }
33
+ ( $ty: ty, $uty: ty, $fraction_bits: expr, $exponent_bias: expr,
34
+ $next_u: ident) => {
35
+ impl IntoFloat for $uty {
36
+ type F = $ty;
37
+ #[ inline( always) ]
38
+ fn into_float_with_exponent( self , exponent: i32 ) -> $ty {
39
+ // The exponent is encoded using an offset-binary representation
40
+ let exponent_bits =
41
+ ( ( $exponent_bias + exponent) as $uty) << $fraction_bits;
42
+ unsafe { mem:: transmute( self | exponent_bits) }
114
43
}
115
- impl Distribution <$ty> for Closed01 {
116
- #[ inline]
117
- fn sample<R : Rng + ?Sized >( & self , rng: & mut R ) -> $ty {
118
- // rescale so that 1.0 - epsilon becomes 1.0
119
- // precisely.
120
- let x: $ty = rng. gen ( ) ;
121
- x * SCALE / ( SCALE - 1.0 )
122
- }
44
+ }
45
+
46
+ impl Distribution <$ty> for Uniform {
47
+ /// Generate a floating point number in the open interval `(0, 1)`
48
+ /// (not including either endpoint) with a uniform distribution.
49
+ ///
50
+ /// # Example
51
+ /// ```rust
52
+ /// use rand::{weak_rng, Rng};
53
+ /// use rand::distributions::Uniform;
54
+ ///
55
+ /// let val: f32 = weak_rng().sample(Uniform);
56
+ /// println!("f32 from (0,1): {}", val);
57
+ /// ```
58
+ fn sample<R : Rng + ?Sized >( & self , rng: & mut R ) -> $ty {
59
+ const EPSILON : $ty = 1.0 / ( 1u64 << $fraction_bits) as $ty;
60
+ let float_size = mem:: size_of:: <$ty>( ) * 8 ;
61
+
62
+ let value = rng. $next_u( ) ;
63
+ let fraction = value >> ( float_size - $fraction_bits) ;
64
+ fraction. into_float_with_exponent( 0 ) - ( 1.0 - EPSILON / 2.0 )
123
65
}
124
66
}
125
67
}
126
68
}
127
- float_impls ! { f64_rand_impls , f64 , 52 }
128
- float_impls ! { f32_rand_impls , f32 , 23 }
69
+ float_impls ! { f32 , u32 , 23 , 127 , next_u32 }
70
+ float_impls ! { f64 , u64 , 52 , 1023 , next_u64 }
129
71
130
72
131
73
#[ cfg( test) ]
132
74
mod tests {
133
75
use Rng ;
134
76
use mock:: StepRng ;
135
- use distributions:: { Open01 , Closed01 } ;
136
77
137
78
const EPSILON32 : f32 = :: core:: f32:: EPSILON ;
138
79
const EPSILON64 : f64 = :: core:: f64:: EPSILON ;
139
80
140
81
#[ test]
141
82
fn floating_point_edge_cases ( ) {
142
83
let mut zeros = StepRng :: new ( 0 , 0 ) ;
143
- assert_eq ! ( zeros. gen :: <f32 >( ) , 0.0 ) ;
144
- assert_eq ! ( zeros. gen :: <f64 >( ) , 0.0 ) ;
145
-
146
- let mut one = StepRng :: new ( 1 , 0 ) ;
147
- assert_eq ! ( one. gen :: <f32 >( ) , EPSILON32 ) ;
148
- assert_eq ! ( one. gen :: <f64 >( ) , EPSILON64 ) ;
149
-
150
- let mut max = StepRng :: new ( !0 , 0 ) ;
151
- assert_eq ! ( max. gen :: <f32 >( ) , 1.0 - EPSILON32 ) ;
152
- assert_eq ! ( max. gen :: <f64 >( ) , 1.0 - EPSILON64 ) ;
153
- }
84
+ assert_eq ! ( zeros. gen :: <f32 >( ) , 0.0 + EPSILON32 / 2.0 ) ;
85
+ assert_eq ! ( zeros. gen :: <f64 >( ) , 0.0 + EPSILON64 / 2.0 ) ;
154
86
155
- #[ test]
156
- fn fp_closed_edge_cases ( ) {
157
- let mut zeros = StepRng :: new ( 0 , 0 ) ;
158
- assert_eq ! ( zeros. sample:: <f32 , _>( Closed01 ) , 0.0 ) ;
159
- assert_eq ! ( zeros. sample:: <f64 , _>( Closed01 ) , 0.0 ) ;
160
-
161
- let mut one = StepRng :: new ( 1 , 0 ) ;
162
- let one32 = one. sample :: < f32 , _ > ( Closed01 ) ;
163
- let one64 = one. sample :: < f64 , _ > ( Closed01 ) ;
164
- assert ! ( EPSILON32 < one32 && one32 < EPSILON32 * 1.01 ) ;
165
- assert ! ( EPSILON64 < one64 && one64 < EPSILON64 * 1.01 ) ;
166
-
167
- let mut max = StepRng :: new ( !0 , 0 ) ;
168
- assert_eq ! ( max. sample:: <f32 , _>( Closed01 ) , 1.0 ) ;
169
- assert_eq ! ( max. sample:: <f64 , _>( Closed01 ) , 1.0 ) ;
170
- }
171
-
172
- #[ test]
173
- fn fp_open_edge_cases ( ) {
174
- let mut zeros = StepRng :: new ( 0 , 0 ) ;
175
- assert_eq ! ( zeros. sample:: <f32 , _>( Open01 ) , 0.0 + EPSILON32 / 2.0 ) ;
176
- assert_eq ! ( zeros. sample:: <f64 , _>( Open01 ) , 0.0 + EPSILON64 / 2.0 ) ;
177
-
178
- let mut one = StepRng :: new ( 1 , 0 ) ;
179
- let one32 = one. sample :: < f32 , _ > ( Open01 ) ;
180
- let one64 = one. sample :: < f64 , _ > ( Open01 ) ;
87
+ let mut one = StepRng :: new ( 1 << 9 , 0 ) ;
88
+ let one32 = one. gen :: < f32 > ( ) ;
181
89
assert ! ( EPSILON32 < one32 && one32 < EPSILON32 * 2.0 ) ;
182
- assert ! ( EPSILON64 < one64 && one64 < EPSILON64 * 2.0 ) ;
183
-
184
- let mut max = StepRng :: new ( !0 , 0 ) ;
185
- assert_eq ! ( max. sample:: <f32 , _>( Open01 ) , 1.0 - EPSILON32 / 2.0 ) ;
186
- assert_eq ! ( max. sample:: <f64 , _>( Open01 ) , 1.0 - EPSILON64 / 2.0 ) ;
187
- }
188
90
189
- #[ test]
190
- fn rand_open ( ) {
191
- // this is unlikely to catch an incorrect implementation that
192
- // generates exactly 0 or 1, but it keeps it sane.
193
- let mut rng = :: test:: rng ( 510 ) ;
194
- for _ in 0 ..1_000 {
195
- // strict inequalities
196
- let f: f64 = rng. sample ( Open01 ) ;
197
- assert ! ( 0.0 < f && f < 1.0 ) ;
198
-
199
- let f: f32 = rng. sample ( Open01 ) ;
200
- assert ! ( 0.0 < f && f < 1.0 ) ;
201
- }
202
- }
203
-
204
- #[ test]
205
- fn rand_closed ( ) {
206
- let mut rng = :: test:: rng ( 511 ) ;
207
- for _ in 0 ..1_000 {
208
- // strict inequalities
209
- let f: f64 = rng. sample ( Closed01 ) ;
210
- assert ! ( 0.0 <= f && f <= 1.0 ) ;
91
+ let mut one = StepRng :: new ( 1 << 12 , 0 ) ;
92
+ let one64 = one. gen :: < f64 > ( ) ;
93
+ assert ! ( EPSILON64 < one64 && one64 < EPSILON64 * 2.0 ) ;
211
94
212
- let f : f32 = rng . sample ( Closed01 ) ;
213
- assert ! ( 0.0 <= f && f <= 1 .0) ;
214
- }
95
+ let mut max = StepRng :: new ( ! 0 , 0 ) ;
96
+ assert_eq ! ( max . gen :: < f32 > ( ) , 1.0 - EPSILON32 / 2 .0) ;
97
+ assert_eq ! ( max . gen :: < f64 > ( ) , 1.0 - EPSILON64 / 2.0 ) ;
215
98
}
216
99
}
0 commit comments