@@ -24,7 +24,7 @@ type MersenneTwister <: AbstractRNG
24
24
MersenneTwister (seed= 0 ) = MersenneTwister (make_seed (seed))
25
25
end
26
26
27
- function srand (r:: MersenneTwister , seed)
27
+ function srand (r:: MersenneTwister , seed)
28
28
r. seed = seed
29
29
dsfmt_init_gen_rand (r. state, seed)
30
30
return r
@@ -146,77 +146,87 @@ rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try
146
146
rand {T<:Number} (:: Type{T} , dims:: Int... ) = rand (T, dims)
147
147
148
148
149
- # # Generate random integer within a range 1:n
149
+ # # Generate random integer within a range
150
150
151
- # maxmultiple: precondition: k>0
152
- # maximum multiple of k <= 2^numbits(k), decremented by one
153
- maxmultiple (k:: Uint32 ) = itrunc (Uint32, div (0x0000000100000000 , k)* k - 1 )
154
- maxmultiple (k:: Uint64 ) = itrunc (Uint64, div (0x00000000000000010000000000000000 , k)* k - 1 )
155
- # maximum multiple of k <= typemax(Uint128), decremented by one
156
- maxmultiple (k:: Uint128 ) = div (typemax (Uint128), k)* k - 1
157
- # 32 or 64 bits version depending on size
158
- maxmultiplemix (k:: Uint64 ) = k >> 32 == 0 ? uint64 (maxmultiple (itrunc (Uint32, k))) : maxmultiple (k)
151
+ # remainder function according to Knuth, where rem_knuth(a, 0) = a
152
+ rem_knuth (a:: Uint , b:: Uint ) = a % (b + (b == 0 )) + a * (b == 0 )
153
+ rem_knuth {T<:Unsigned} (a:: T , b:: T ) = b != 0 ? a % b : a
159
154
160
- immutable RandIntGen{T<: Integer , U<: Union(Uint32, Uint64, Uint128) }
161
- k:: U # range length
162
- u:: U # rejection threshold
155
+ # maximum multiple of k <= 2^bits(T) decremented by one,
156
+ # that is 0xFFFFFFFF if k = typemax(T) - typemin(T) with intentional underflow
157
+ maxmultiple (k:: Uint32 ) = itrunc (Uint32, div (0x0000000100000000 ,k + (k == 0 ))* k - 1 )
158
+ maxmultiple (k:: Uint64 ) = itrunc (Uint64, div (0x00000000000000010000000000000000 , k + (k == 0 ))* k - 1 )
159
+ # maximum multiple of k within 1:typemax(Uint128)
160
+ maxmultiple (k:: Uint128 ) = div (typemax (Uint128), k + (k == 0 ))* k - 1
161
+ # maximum multiple of k within 1:2^32 or 1:2^64, depending on size
162
+ maxmultiplemix (k:: Uint64 ) = itrunc (Uint64, div ((k >> 32 != 0 )* 0x0000000000000000FFFFFFFF00000000 + 0x0000000100000000 , k + (k == 0 ))* k - 1 )
163
163
164
- function RandIntGen (k:: U )
165
- k < 1 && error (" No integers in the range [1, $k ]. Did you try to pick an element from a empty collection?" )
166
- new (k, isa (k, Uint64) ? maxmultiplemix (k) : maxmultiple (k))
167
- end
164
+ immutable RandIntGen{T<: Integer , U<: Unsigned }
165
+ k:: U # range length (or zero for full range)
166
+ u:: U # rejection threshold
168
167
end
168
+ # generators with 32, 128 bits entropy
169
+ RandIntGen {U<:Union(Uint32, Uint128)} (T:: Type , k:: U ) = RandIntGen {T, U} (k, maxmultiple (k))
170
+ # mixed 32/64 bits entropy generator
171
+ RandIntGen (T:: Type , k:: Uint64 ) = RandIntGen {T,Uint64} (k, maxmultiplemix (k))
172
+
169
173
170
174
# generator API
171
175
# randintgen(k) returns a helper object for generating random integers in the range 1:k
172
- randintgen {T<:Unsigned} (k:: T ) = RandIntGen { T, T} ( k)
176
+ randintgen {T<:Unsigned} (k:: T ) = k < 1 ? error ( " range must be non-empty " ) : RandIntGen ( T, k)
173
177
# specialized versions
174
178
for (T, U) in [(Uint8, Uint32), (Uint16, Uint32),
175
179
(Int8, Uint32), (Int16, Uint32), (Int32, Uint32), (Int64, Uint64), (Int128, Uint128),
176
180
(Bool, Uint32), (Char, Uint32)]
177
181
178
- @eval randintgen (k:: $T ) = RandIntGen {$T,$U} (convert ($ U, k))
179
- end
180
-
181
- @inline function rand_lessthan {U} (u:: U )
182
- while true
183
- x = rand (U)
184
- x <= u && return x
185
- end
182
+ @eval randintgen (k:: $T ) = k< 1 ? error (" range must be non-empty" ) : RandIntGen ($ T, convert ($ U, k)) # overflow ok
186
183
end
187
184
188
185
# this function uses 32 bit entropy for small ranges of length <= typemax(Uint32) + 1
189
186
# RandIntGen is responsible for providing the right value of k
190
- function rand {T} (g:: RandIntGen{T,Uint64} )
191
- x = (g. k - 1 ) >> 32 == 0 ?
192
- uint64 (rand_lessthan (itrunc (Uint32, g. u))) :
193
- rand_lessthan (g. u)
194
- return reinterpret (T, one (Uint64) + x % g. k)
187
+ function rand {T<:Union(Uint64, Int64)} (g:: RandIntGen{T,Uint64} )
188
+ local x:: Uint64
189
+ if (g. k - 1 ) >> 32 == 0
190
+ x = rand (Uint32)
191
+ while x > g. u
192
+ x = rand (Uint32)
193
+ end
194
+ else
195
+ x = rand (Uint64)
196
+ while x > g. u
197
+ x = rand (Uint64)
198
+ end
199
+ end
200
+ return reinterpret (T, one (Uint64) + rem_knuth (x, g. k))
195
201
end
196
202
197
- rand {T,U } (g:: RandIntGen{T,U} ) = itrunc (T, one (U) + rand_lessthan (g . u) % g . k )
198
-
199
-
200
- # # Randomly draw a sample from an AbstractArray r
201
- # (e.g. r is a range 0:2:8 or a vector [2, 3, 5, 7])
202
-
203
- rand (r :: AbstractArray ) = r[ rand ( randintgen ( length (r)))]
203
+ function rand {T<:Integer, U<:Unsigned } (g:: RandIntGen{T,U} )
204
+ x = rand (U)
205
+ while x > g . u
206
+ x = rand (U)
207
+ end
208
+ itrunc (T, one (U) + rem_knuth (x, g . k))
209
+ end
204
210
205
- # Arrays of random integers
211
+ rand {T} (r :: Range{T} ) = r[ rand ( randintgen ( length (r)))]
206
212
207
- rand! (r:: Range , A:: AbstractArray ) = rand! (r, A, ())
213
+ function rand! (g:: RandIntGen , A:: AbstractArray )
214
+ for i = 1 : length (A)
215
+ @inbounds A[i] = rand (g)
216
+ end
217
+ return A
218
+ end
208
219
209
- # TODO : this more general version is "disabled" until #8246 is resolved
210
- function rand! (r:: AbstractArray , A:: AbstractArray , :: () )
220
+ function rand! (r:: Range , A:: AbstractArray )
211
221
g = randintgen (length (r))
212
222
for i = 1 : length (A)
213
223
@inbounds A[i] = r[rand (g)]
214
224
end
215
225
return A
216
226
end
217
227
218
- rand {T} (r:: AbstractArray {T} , dims:: Dims ) = rand! (r, Array (T, dims), ( ))
219
- rand (r:: AbstractArray , dims:: Int... ) = rand (r, dims)
228
+ rand {T} (r:: Range {T} , dims:: Dims ) = rand! (r, Array (T, dims))
229
+ rand (r:: Range , dims:: Int... ) = rand (r, dims)
220
230
221
231
222
232
# # random Bools
0 commit comments