Merge pull request #8854 from JuliaLang/rf/rand-rng

RFC: more methods working with a non-global RNG

base/bitarray.jl

@@ -331,10 +331,10 @@ bitpack{T,N}(A::AbstractArray{T,N}) = convert(BitArray{N}, A)
 
 ## Random ##
 
-function bitarray_rand_fill!(B::BitArray)
+function bitarray_rand_fill!(rng, B::BitArray) # rng is an AbstractRNG
     length(B) == 0 && return B
     Bc = B.chunks
-    rand!(Bc)
+    rand!(rng, Bc)
     Bc[end] &= @_msk_end length(B)
     return B
 end

base/random.jl

@@ -13,17 +13,13 @@ abstract AbstractRNG
 
 type MersenneTwister <: AbstractRNG
     state::DSFMT_state
-    seed::Union(Uint32,Vector{Uint32})
     vals::Vector{Float64}
     idx::Int
+    seed::Vector{Uint32}
 
-    function MersenneTwister(seed::Vector{Uint32})
-        state = DSFMT_state()
-        dsfmt_init_by_array(state, seed)
-        return new(state, seed, Array(Float64, dsfmt_get_min_array_size()), dsfmt_get_min_array_size())
-    end
-
-    MersenneTwister(seed=0) = MersenneTwister(make_seed(seed))
+    MersenneTwister(seed) = srand(new(DSFMT_state(), Array(Float64, dsfmt_get_min_array_size())),
+                                  seed)
+    MersenneTwister() = MersenneTwister(0)
 end
 
 ## Low level API for MersenneTwister
@@ -47,50 +43,43 @@ end
 @inline rand_ui32(r::MersenneTwister) = reinterpret(Uint64, rand_close1_open2(r)) % Uint32
 
 
-function srand(r::MersenneTwister, seed)
+function srand(r::MersenneTwister, seed::Vector{Uint32})
     r.seed = seed
-    dsfmt_init_gen_rand(r.state, seed)
+    dsfmt_init_by_array(r.state, r.seed)
     r.idx = length(r.vals)
     return r
 end
 
 ## initialization
 
-function srand()
+__init__() = srand()
+
+## make_seed()
+# make_seed methods produce values of type Array{Uint32}, suitable for MersenneTwister seeding
+
+function make_seed()
 
 @unix_only begin
     try
-        srand("/dev/urandom")
+        return make_seed("/dev/urandom", 4)
     catch
         println(STDERR, "Entropy pool not available to seed RNG; using ad-hoc entropy sources.")
         seed = reinterpret(Uint64, time())
         seed = hash(seed, uint64(getpid()))
         try
         seed = hash(seed, parseint(Uint64, readall(`ifconfig` |> `sha1sum`)[1:40], 16))
         end
-        srand(seed)
+        return make_seed(seed)
     end
 end
 
 @windows_only begin
     a = zeros(Uint32, 2)
     win32_SystemFunction036!(a)
-    srand(a)
+    return a
 end
 end
 
-__init__() = srand()
-
-## srand()
-
-function srand(seed::Vector{Uint32})
-    GLOBAL_RNG.seed = seed
-    dsfmt_init_by_array(GLOBAL_RNG.state, seed)
-    GLOBAL_RNG.idx = length(GLOBAL_RNG.vals)
-    return GLOBAL_RNG
-end
-srand(n::Integer) = srand(make_seed(n))
-
 function make_seed(n::Integer)
     n < 0 && throw(DomainError())
     seed = Uint32[]
@@ -103,63 +92,76 @@ function make_seed(n::Integer)
     end
 end
 
-function srand(filename::String, n::Integer)
+function make_seed(filename::String, n::Integer)
     open(filename) do io
         a = Array(Uint32, int(n))
         read!(io, a)
-        srand(a)
+        a
     end
 end
-srand(filename::String) = srand(filename, 4)
+
+## srand()
+
+srand(r::MersenneTwister) = srand(r, make_seed())
+srand(r::MersenneTwister, n::Integer) = srand(r, make_seed(n))
+srand(r::MersenneTwister, filename::String, n::Integer=4) = srand(r, make_seed(filename, n))
+
+srand() = srand(GLOBAL_RNG)
+srand(seed::Union(Integer, Vector{Uint32})) = srand(GLOBAL_RNG, seed)
+srand(filename::String, n::Integer=4) = srand(GLOBAL_RNG, filename, n)
 
 ## Global RNG
 
 const GLOBAL_RNG = MersenneTwister()
 globalRNG() = GLOBAL_RNG
 
-## random floating point values
+# rand: a non-specified RNG defaults to GLOBAL_RNG
 
-rand(r::MersenneTwister=GLOBAL_RNG) = rand_close_open(r)
+rand() = rand(GLOBAL_RNG)
+rand(T::Type) = rand(GLOBAL_RNG, T)
+rand(::()) = rand(GLOBAL_RNG, ()) # needed to resolve ambiguity
+rand(dims::Dims) = rand(GLOBAL_RNG, dims)
+rand(dims::Int...) = rand(dims)
+rand(T::Type, dims::Dims) = rand(GLOBAL_RNG, T, dims)
+rand(T::Type, d1::Int, dims::Int...) = rand(T, tuple(d1, dims...))
+rand!(A::AbstractArray) = rand!(GLOBAL_RNG, A)
 
-rand(::Type{Float64}) = rand()
+## random floating point values
 
-rand(::Type{Float32}) = float32(rand())
-rand(::Type{Float16}) = float16(rand())
+rand(r::AbstractRNG) = rand(r, Float64)
 
-rand{T<:Real}(::Type{Complex{T}}) = complex(rand(T),rand(T))
+# MersenneTwister
+rand(r::MersenneTwister, ::Type{Float64}) = rand_close_open(r)
+rand{T<:Union(Float16, Float32)}(r::MersenneTwister, ::Type{T}) = convert(T, rand(r, Float64))
 
-## random integers
+## random integers (MersenneTwister)
 
-rand(::Type{Uint8})   = rand(Uint32) % Uint8
-rand(::Type{Uint16})  = rand(Uint32) % Uint16
-rand(::Type{Uint32})  = rand_ui32(GLOBAL_RNG)
-rand(::Type{Uint64})  = uint64(rand(Uint32)) <<32 | rand(Uint32)
-rand(::Type{Uint128}) = uint128(rand(Uint64))<<64 | rand(Uint64)
+rand(r::MersenneTwister, ::Type{Uint8})   = rand(r, Uint32) % Uint8
+rand(r::MersenneTwister, ::Type{Uint16})  = rand(r, Uint32) % Uint16
+rand(r::MersenneTwister, ::Type{Uint32})  = rand_ui32(r)
+rand(r::MersenneTwister, ::Type{Uint64})  = uint64(rand(r, Uint32)) <<32 | rand(r, Uint32)
+rand(r::MersenneTwister, ::Type{Uint128}) = uint128(rand(r, Uint64))<<64 | rand(r, Uint64)
 
-rand(::Type{Int8})    = rand(Uint32) % Int8
-rand(::Type{Int16})   = rand(Uint32) % Int16
-rand(::Type{Int32})   = reinterpret(Int32,rand(Uint32))
-rand(::Type{Int64})   = reinterpret(Int64,rand(Uint64))
-rand(::Type{Int128})  = reinterpret(Int128,rand(Uint128))
+rand(r::MersenneTwister, ::Type{Int8})    = rand(r, Uint32) % Int8
+rand(r::MersenneTwister, ::Type{Int16})   = rand(r, Uint32) % Int16
+rand(r::MersenneTwister, ::Type{Int32})   = reinterpret(Int32,  rand(r, Uint32))
+rand(r::MersenneTwister, ::Type{Int64})   = reinterpret(Int64,  rand(r, Uint64))
+rand(r::MersenneTwister, ::Type{Int128})  = reinterpret(Int128, rand(r, Uint128))
 
-# Arrays of random numbers
+## random complex values
 
-rand(::Type{Float64}, dims::Dims) = rand!(Array(Float64, dims))
-rand(::Type{Float64}, dims::Int...) = rand(Float64, dims)
+rand{T<:Real}(r::AbstractRNG, ::Type{Complex{T}}) = complex(rand(r, T), rand(r, T))
 
-rand(dims::Dims) = rand(Float64, dims)
-rand(dims::Int...) = rand(Float64, dims)
+## Arrays of random numbers
 
-rand(r::AbstractRNG) = rand(r, Float64)
-rand(r::AbstractRNG, dims::Dims) = rand!(r, Array(Float64, dims))
+rand(r::AbstractRNG, dims::Dims) = rand(r, Float64, dims)
 rand(r::AbstractRNG, dims::Int...) = rand(r, dims)
 
-function rand!{T}(A::Array{T})
-    for i = 1:length(A)
-        A[i] = rand(T)
-    end
-    A
-end
+rand(r::AbstractRNG, T::Type, dims::Dims) = rand!(r, Array(T, dims))
+rand(r::AbstractRNG, T::Type, d1::Int, dims::Int...) = rand(r, T, tuple(d1, dims...))
+# note: the above method would trigger an ambiguity warning if d1 was not separated out:
+# rand(r, ()) would match both this method and rand(r, dims::Dims)
+# moreover, a call like rand(r, NotImplementedType()) would be an infinite loop
 
 function rand!{T}(r::AbstractRNG, A::AbstractArray{T})
     for i = 1:length(A)
@@ -168,6 +170,8 @@ function rand!{T}(r::AbstractRNG, A::AbstractArray{T})
     A
 end
 
+# MersenneTwister
+
 function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float64})
     n = length(A)
     # what follows is equivalent to this simple loop but more efficient:
@@ -203,14 +207,6 @@ function rand!(r::MersenneTwister, A::Array{Float64})
     A
 end
 
-rand!(A::AbstractArray{Float64}) = rand!(GLOBAL_RNG, A)
-rand!(A::Array{Float64}) = rand!(GLOBAL_RNG, A)
-
-rand(T::Type, dims::Dims) = rand!(Array(T, dims))
-rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try a more specific type")
-rand{T<:Number}(::Type{T}, dims::Int...) = rand(T, dims)
-
-
 ## Generate random integer within a range
 
 # remainder function according to Knuth, where rem_knuth(a, 0) = a
@@ -297,16 +293,20 @@ end
 rand{T}(r::Range{T}, dims::Dims) = rand!(r, Array(T, dims))
 rand(r::Range, dims::Int...) = rand(r, dims)
 
+## random BitArrays (AbstractRNG)
 
-## random Bools
+rand!(r::AbstractRNG, B::BitArray) = Base.bitarray_rand_fill!(r, B)
 
-rand!(B::BitArray) = Base.bitarray_rand_fill!(B)
+randbool(r::AbstractRNG, dims::Dims)   = rand!(r, BitArray(dims))
+randbool(r::AbstractRNG, dims::Int...) = rand!(r, BitArray(dims))
 
-randbool(dims::Dims) = rand!(BitArray(dims))
+randbool(dims::Dims)   = rand!(BitArray(dims))
 randbool(dims::Int...) = rand!(BitArray(dims))
 
-randbool() = ((rand(Uint32) & 1) == 1)
-rand(::Type{Bool}) = randbool()
+randbool(r::MersenneTwister=GLOBAL_RNG) = ((rand(r, Uint32) & 1) == 1)
+
+rand(r::MersenneTwister, ::Type{Bool}) = randbool(r)
+
 
 ## randn() - Normally distributed random numbers using Ziggurat algorithm
 

doc/stdlib/base.rst

@@ -3931,43 +3931,35 @@ Random Numbers
 
 Random number generation in Julia uses the `Mersenne Twister library <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT>`_. Julia has a global RNG, which is used by default. Multiple RNGs can be plugged in using the ``AbstractRNG`` object, which can then be used to have multiple streams of random numbers. Currently, only ``MersenneTwister`` is supported.
 
+Most functions related to random generation accept an optional ``AbstractRNG`` as the first argument, ``rng`` , which defaults to the global one if not provided. Morever, some of them accept optionally dimension specifications ``dims...`` (which can be given as a tuple) to generate arrays of random values.
+
+A ``MersenneTwister`` RNG can generate random numbers of the following types: ``Float16, Float32, Float64, Bool, Int16, Uint16, Int32, Uint32, Int64, Uint64, Int128, Uint128`` (or complex numbers or arrays of those types). Random floating point numbers are generated uniformly in [0,1).
+
 .. function:: srand([rng], [seed])
 
-   Reseed the random number generator. If a ``seed`` is provided, the RNG will give a reproducible sequence of numbers, otherwise Julia will get entropy from the system. The ``seed`` may be an unsigned integer, a vector of unsigned integers or a filename, in which case the seed is read from a file. If the argument ``rng`` is not provided, the default global RNG is seeded.
+   Reseed the random number generator. If a ``seed`` is provided, the RNG will give a reproducible sequence of numbers, otherwise Julia will get entropy from the system. The ``seed`` may be a non-negative integer, a vector of ``Uint32`` integers or a filename, in which case the seed is read from a file.
 
 .. function:: MersenneTwister([seed])
 
    Create a ``MersenneTwister`` RNG object. Different RNG objects can have their own seeds, which may be useful for generating different streams of random numbers.
 
-.. function:: rand() -> Float64
+.. function:: rand([rng], [t::Type], [dims...])
 
-   Generate a ``Float64`` random number uniformly in [0,1)
+   Generate a random value or an array of random values of the given type, ``t``, which defaults to ``Float64``.
 
 .. function:: rand!([rng], A)
 
-   Populate the array A with random number generated from the specified RNG.
-
-.. function:: rand(rng::AbstractRNG, [dims...])
-
-   Generate a random ``Float64`` number or array of the size specified by dims, using the specified RNG object. Currently, ``MersenneTwister`` is the only available Random Number Generator (RNG), which may be seeded using srand.
-
-.. function:: rand(dims or [dims...])
-
-   Generate a random ``Float64`` array of the size specified by dims
-
-.. function:: rand(t::Type, [dims...])
-
-   Generate a random number or array of random numbes of the given type.
+   Populate the array A with random values.
 
 .. function:: rand(r, [dims...])
 
    Pick a random element or array of random elements from range ``r`` (for example, ``1:n`` or ``0:2:10``).
 
-.. function:: randbool([dims...])
+.. function:: randbool([rng], [dims...])
 
-   Generate a random boolean value. Optionally, generate an array of random boolean values.
+   Generate a random boolean value. Optionally, generate a ``BitArray`` of random boolean values.
 
-.. function:: randn([rng], dims or [dims...])
+.. function:: randn([rng], [dims...])
 
    Generate a normally-distributed random number with mean 0 and standard deviation 1. Optionally generate an array of normally-distributed random numbers.
 

test/random.jl

@@ -19,12 +19,17 @@ srand(0); rand(); x = rand(384);
 # Try a seed larger than 2^32
 @test rand(MersenneTwister(5294967296)) == 0.3498809918210497
 
-# Test array filling, Issue #7643
+# Test array filling, Issues #7643, #8360
 @test rand(MersenneTwister(0), 1) == [0.8236475079774124]
 A = zeros(2, 2)
 rand!(MersenneTwister(0), A)
 @test A == [0.8236475079774124  0.16456579813368521;
             0.9103565379264364  0.17732884646626457]
+@test rand(MersenneTwister(0), Int64, 1) == [172014471070449746]
+A = zeros(Int64, 2, 2)
+rand!(MersenneTwister(0), A)
+@test A == [  172014471070449746  -193283627354378518;
+            -4679130500738884555 -9008350441255501549]
 
 # randn
 @test randn(MersenneTwister(42)) == -0.5560268761463861