rand of AbstractArray

base/random.jl

@@ -1,6 +1,6 @@
 module Random
 
-using Base.dSFMT
+using Base: dSFMT, IntTypes
 
 export srand,
        rand, rand!,
@@ -146,91 +146,88 @@ rand{T<:Number}(::Type{T}) = error("no random number generator for type $T; try
 rand{T<:Number}(::Type{T}, dims::Int...) = rand(T, dims)
 
 
-## Generate random integer within a range
+## Generate random integer within a range 0:k-1
 
-# remainder function according to Knuth, where rem_knuth(a, 0) = a
-rem_knuth(a::Uint, b::Uint) = a % (b + (b == 0)) + a * (b == 0)
-rem_knuth{T<:Unsigned}(a::T, b::T) = b != 0 ? a % b : a
+# maxmultiple: precondition: k>0
+# maximum multiple of k <= 2^numbits(k), decremented by one
+maxmultiple(k::Uint32) = itrunc(Uint32, div(0x0000000100000000, k)*k - 1)
+maxmultiple(k::Uint64) = itrunc(Uint64, div(0x00000000000000010000000000000000, k)*k - 1)
+# maximum multiple of k <= typemax(Uint128), decremented by one
+maxmultiple(k::Uint128) = div(typemax(Uint128), k)*k - 1
+# 32 or 64 bits version depending on size
+maxmultiplemix(k::Uint64) = k >> 32 == 0 ? uint64(maxmultiple(itrunc(Uint32, k))) : maxmultiple(k)
 
-# maximum multiple of k <= 2^bits(T) decremented by one,
-# that is 0xFFFFFFFF if k = typemax(T) - typemin(T) with intentional underflow
-maxmultiple(k::Uint32) = itrunc(Uint32, div(0x0000000100000000,k + (k == 0))*k - 1)
-maxmultiple(k::Uint64) = itrunc(Uint64, div(0x00000000000000010000000000000000, k + (k == 0))*k - 1)
-# maximum multiple of k within 1:typemax(Uint128)
-maxmultiple(k::Uint128) = div(typemax(Uint128), k + (k == 0))*k - 1
-# maximum multiple of k within 1:2^32 or 1:2^64, depending on size
-maxmultiplemix(k::Uint64) = itrunc(Uint64, div((k >> 32 != 0)*0x0000000000000000FFFFFFFF00000000 + 0x0000000100000000, k + (k == 0))*k - 1)
-
-immutable RandIntGen{T<:Integer, U<:Unsigned}
-    a::T   # first element of the range
+immutable RandIntGen{U<:Union(Uint32, Uint64, Uint128)}
     k::U   # range length or zero for full range
     u::U   # rejection threshold
-end
-# generators with 32, 128 bits entropy
-RandIntGen{T, U<:Union(Uint32, Uint128)}(a::T, k::U) = RandIntGen{T, U}(a, k, maxmultiple(k))
-# mixed 32/64 bits entropy generator
-RandIntGen{T}(a::T, k::Uint64) = RandIntGen{T,Uint64}(a, k, maxmultiplemix(k))
-
 
-# generator for ranges
-RandIntGen{T<:Unsigned}(r::UnitRange{T}) = isempty(r) ? error("range must be non-empty") : RandIntGen(first(r), last(r) - first(r) + one(T))
+    function RandIntGen(k::U)
+        u = k == zero(U) ? typemax(U) :
+            U === Uint64 ? maxmultiplemix(k) :
+                           maxmultiple(k)
+        new(k, u)
+    end
+end
 
+# generator API
+# randintgen(k) returns a helper object for generating random integers in the range 0:k
+randintgen{U<:Unsigned}(k::U) = RandIntGen{U}(one(U)+k)
 # specialized versions
-for (T, U) in [(Uint8, Uint32), (Uint16, Uint32),
-               (Int8, Uint32), (Int16, Uint32), (Int32, Uint32), (Int64, Uint64), (Int128, Uint128),
-               (Bool, Uint32), (Char, Uint32)]
+randintgen(k::Uint8)  = randintgen(convert(Uint32, k))
+randintgen(k::Uint16) = randintgen(convert(Uint32, k))
 
-    @eval RandIntGen(r::UnitRange{$T}) = isempty(r) ? error("range must be non-empty") : RandIntGen(first(r), convert($U, unsigned(last(r) - first(r)) + one($U))) # overflow ok
+@inline function rand_lessthan{U}(u::U)
+    while true
+        x = rand(U)
+        x <= u && return x
+    end
 end
 
 # this function uses 32 bit entropy for small ranges of length <= typemax(Uint32) + 1
 # RandIntGen is responsible for providing the right value of k
-function rand{T<:Union(Uint64, Int64)}(g::RandIntGen{T,Uint64})
-    local x::Uint64
-    if (g.k - 1) >> 32 == 0
-        x = rand(Uint32)
-        while x > g.u
-            x = rand(Uint32)
-        end
-    else
-        x = rand(Uint64)
-        while x > g.u
-            x = rand(Uint64)
-        end
-    end
-    return reinterpret(T, reinterpret(Uint64, g.a) + rem_knuth(x, g.k))
+function rand(g::RandIntGen{Uint64})
+    g.k == zero(Uint64) && return rand(Uint64)
+    x = (g.k - 1) >> 32 == 0 ?
+        uint64(rand_lessthan(itrunc(Uint32, g.u))) :
+        rand_lessthan(g.u)
+    return x % g.k
 end
 
-function rand{T<:Integer, U<:Unsigned}(g::RandIntGen{T,U})
-    x = rand(U)
-    while x > g.u
-        x = rand(U)
-    end
-    itrunc(T, unsigned(g.a) + rem_knuth(x, g.k))
-end
+rand{U}(g::RandIntGen{U}) = g.k == zero(U) ? rand(U) : rand_lessthan(g.u) % g.k
 
-rand{T<:Union(Signed,Unsigned,Bool,Char)}(r::UnitRange{T}) = rand(RandIntGen(r))
-rand{T}(r::Range{T}) = r[rand(1:(length(r)))]
 
-function rand!(g::RandIntGen, A::AbstractArray)
-    for i = 1 : length(A)
-        @inbounds A[i] = rand(g)
-    end
-    return A
-end
+## Randomly draw a sample from an AbstractArray r
+# (e.g. r is a range 0:2:8 or a vector [2, 3, 5, 7])
+
+check_nonempty(r::AbstractArray) = (isempty(r) && error(ArgumentError("cannot randomly pick an element from an empty collection")); r)
+
+# special code which can handle "full ranges" (a workaround to error in e.g. `length(typemin(Int):typemax(Int))`)
+# getindex_0/endof_0 work like getindex/endof respectively, but with indexes starting at 0
+getindex_0(t::AbstractArray, i::Real) = @inbounds return t[i+1]
+getindex_0{T<:Union(IntTypes...)}(r::Range{T}, i::Integer) = itrunc(T, unsigned(first(r)) +  unsigned(i)*unsigned(step(r)))
+endof_0(a::AbstractArray) = unsigned(length(a) - 1) # precondition: !isempty(a)
+endof_0{T<:Union(IntTypes...),S<:Union(IntTypes...)}(r::StepRange{T,S}) = r.step > 0 ?
+                                                                         div(unsigned(r.stop - r.start), r.step) :
+                                                                         div(unsigned(r.start - r.stop), -r.step)
+endof_0{T<:Union(IntTypes...)}(r::UnitRange{T}) = unsigned(r.stop - r.start)
+
+rand(r::AbstractArray) = getindex_0(r, rand(randintgen(endof_0(check_nonempty(r)))))
+
+# Arrays of random integers
 
-rand!{T<:Union(Signed,Unsigned,Bool,Char)}(r::UnitRange{T}, A::AbstractArray) = rand!(RandIntGen(r), A)
+rand!(r::Range, A::AbstractArray) = rand!(r, A, ())
 
-function rand!(r::Range, A::AbstractArray)
-    g = RandIntGen(1:(length(r)))
+# TODO: this more general version is "disabled" until #8246 is resolved
+function rand!(r::AbstractArray, A::AbstractArray, ::())
+    g = randintgen(endof_0(check_nonempty(r)))
     for i = 1 : length(A)
-        @inbounds A[i] = r[rand(g)]
+        @inbounds A[i] = getindex_0(r, rand(g))
     end
     return A
 end
 
-rand{T}(r::Range{T}, dims::Dims) = rand!(r, Array(T, dims))
-rand(r::Range, dims::Int...) = rand(r, dims)
+rand{T}(r::AbstractArray{T}, dims::Dims) = rand!(r, Array(T, dims), ())
+rand(r::AbstractArray, dims::Int...) = rand(r, dims)
 
 
 ## random Bools

doc/stdlib/base.rst

@@ -3947,7 +3947,7 @@ Random number generation in Julia uses the `Mersenne Twister library <http://www
 
 .. function:: rand!([rng], A)
 
-   Populate the array A with random number generated from the specified RNG.
+   Populate the array A with random numbers generated from the specified RNG.
 
 .. function:: rand(rng::AbstractRNG, [dims...])
 
@@ -3961,9 +3961,13 @@ Random number generation in Julia uses the `Mersenne Twister library <http://www
 
    Generate a random number or array of random numbes of the given type.
 
-.. function:: rand(r, [dims...])
+.. function:: rand(coll, [dims...])
 
-   Pick a random element or array of random elements from range ``r`` (for example, ``1:n`` or ``0:2:10``).
+   Pick a random element or array of random elements from the indexable collection  ``coll`` (for example, ``1:n`` or ``['x','y','z']``).
+
+.. function:: rand!(r, A)
+
+   Populate the array A with random values drawn uniformly from the range ``r``.
 
 .. function:: randbool([dims...])
 

test/random.jl

@@ -26,6 +26,11 @@ rand!(MersenneTwister(0), A)
 @test A == [0.8236475079774124  0.16456579813368521;
             0.9103565379264364  0.17732884646626457]
 
+@test rand(0:3:1000) in 0:3:1000
+coll = Any[2, uint128(128), big(619), "string", 'c']
+@test rand(coll) in coll
+@test issubset(rand(coll, 2, 3), coll)
+
 # randn
 @test randn(MersenneTwister(42)) == -0.5560268761463861
 A = zeros(2, 2)
@@ -54,9 +59,7 @@ if sizeof(Int32) < sizeof(Int)
     r = rand(int32(-1):typemax(Int32))
     @test typeof(r) == Int32
     @test -1 <= r <= typemax(Int32)
-    @test all([div(0x00010000000000000000,k)*k - 1 == Base.Random.RandIntGen(uint64(1:k)).u for k in 13 .+ int64(2).^(32:62)])
-    @test all([div(0x00010000000000000000,k)*k - 1 == Base.Random.RandIntGen(int64(1:k)).u for k in 13 .+ int64(2).^(32:61)])
-
+    @test all([div(0x00010000000000000000,k)*k - 1 == Base.Random.randintgen(uint64(k)-1).u for k in 13 .+ int64(2).^(32:62)])
 end
 
 randn(100000)
@@ -160,8 +163,7 @@ r = uint64(rand(uint32(97:122)))
 srand(seed)
 @test r == rand(uint64(97:122))
 
-@test all([div(0x000100000000,k)*k - 1 == Base.Random.RandIntGen(uint64(1:k)).u for k in 13 .+ int64(2).^(1:30)])
-@test all([div(0x000100000000,k)*k - 1 == Base.Random.RandIntGen(int64(1:k)).u for k in 13 .+ int64(2).^(1:30)])
+@test all([div(0x000100000000,k)*k - 1 == Base.Random.randintgen(uint64(k)-1).u for k in 13 .+ int64(2).^(1:30)])
 
 import Base.Random: uuid4, UUID
 

test/ranges.jl

@@ -375,3 +375,16 @@ end
 
 # issue #8584
 @test (0:1//2:2)[1:2:3] == 0:1//1:1
+
+# endof_0, getindex_0
+using Base.Random: endof_0, getindex_0
+for r in (typemin(Int):typemax(Int), typemin(Int):1:typemax(Int))
+    @test endof_0(r)  == typemax(Uint)
+    @test getindex_0(r, 0) == typemin(Int)
+    @test getindex_0(r, endof_0(r)) == typemax(Int)
+end
+let r = typemax(Int):-1:typemin(Int)
+    @test endof_0(r) == typemax(Uint)
+    @test getindex_0(r, 0) == typemax(Int)
+    @test getindex_0(r, endof_0(r)) == typemin(Int)
+end