implement rand(::Range{BigInt})

* Previously, calling e.g. rand(big(1:4)) caused a stack overflow,
  because rand(mt::MersenneTwister, r::AbstractArray) was calling
  itself recursively. So a a method handling not-implemented types
  has been added.
* A type RandIntGenBigInt similar to RandIntGen (both subtypes of
  RangeGenerator) has been introduced to handle rand on big ranges.
  The generic function to construct such objects is named inrange.
  Note: mpz_import could be replaced by mpz_limbs_{write,finish}
        when GMP version 6 is used.
* BigInt tests from commit bf8c452 were re-added.

Author: rfourquet

base/random.jl

@@ -361,7 +361,9 @@ 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) = (div((k >> 32 != 0)*0x0000000000000000FFFFFFFF00000000 + 0x0000000100000000, k + (k == 0))*k - 1) % UInt64
 
-immutable RandIntGen{T<:Integer, U<:Unsigned}
+abstract RangeGenerator
+
+immutable RandIntGen{T<:Integer, U<:Unsigned} <: RangeGenerator
     a::T   # first element of the range
     k::U   # range length or zero for full range
     u::U   # rejection threshold
@@ -373,16 +375,36 @@ 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))
+inrange{T<:Unsigned}(r::UnitRange{T}) = isempty(r) ? error("range must be non-empty") : RandIntGen(first(r), last(r) - first(r) + one(T))
 
 # specialized versions
 for (T, U) in [(UInt8, UInt32), (UInt16, UInt32),
                (Int8, UInt32), (Int16, UInt32), (Int32, UInt32), (Int64, UInt64), (Int128, UInt128),
                (Bool, UInt32)]
 
-    @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
+    @eval inrange(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
+end
+
+
+# generator for BigInt
+immutable RandIntGenBigInt <: RangeGenerator
+    a::BigInt             # first
+    m::BigInt             # range length - 1
+    limbs::Array{Culong}  # buffer to be copied into generated BigInt's
+    mask::Culong          # applied to the highest limb
 end
 
+function inrange(r::UnitRange{BigInt})
+    m = last(r) - first(r)
+    m < 0 && error("range must be non-empty")
+    nd = ndigits(m, 2)
+    nlimbs, highbits = divrem(nd, 8*sizeof(Culong))
+    highbits > 0 && (nlimbs += 1)
+    mask = highbits == 0 ? ~zero(Culong) : one(Culong)<<highbits - one(Culong)
+    return RandIntGenBigInt(first(r), m, Array(Culong, nlimbs), mask)
+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)}(mt::MersenneTwister, g::RandIntGen{T,UInt64})
@@ -409,23 +431,41 @@ function rand{T<:Integer, U<:Unsigned}(mt::MersenneTwister, g::RandIntGen{T,U})
     (unsigned(g.a) + rem_knuth(x, g.k)) % T
 end
 
-rand{T<:Union(Signed,Unsigned,Bool,Char)}(mt::MersenneTwister, r::UnitRange{T}) = rand(mt, RandIntGen(r))
+function rand(rng::AbstractRNG, g::RandIntGenBigInt)
+    x = BigInt()
+    while true
+        rand!(rng, g.limbs)
+        g.limbs[end] &= g.mask
+        ccall((:__gmpz_import, :libgmp), Void,
+              (Ptr{BigInt}, Csize_t, Cint, Csize_t, Cint, Csize_t, Ptr{Void}),
+              &x, length(g.limbs), -1, sizeof(Culong), 0, 0, &g.limbs)
+        x <= g.m && break
+    end
+    ccall((:__gmpz_add, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &x, &x, &g.a)
+    return x
+end
+
+
+rand{T<:Union(Signed,Unsigned,BigInt,Bool,Char)}(mt::MersenneTwister, r::UnitRange{T}) = rand(mt, inrange(r))
+
+rand{T}(mt::MersenneTwister, r::UnitRange{T}) = error("there is no random generator for ranges of type $T")
+
 
 # Randomly draw a sample from an AbstractArray r
 # (e.g. r is a range 0:2:8 or a vector [2, 3, 5, 7])
 rand(mt::MersenneTwister, r::AbstractArray) = @inbounds return r[rand(mt, 1:length(r))]
 
-function rand!(mt::MersenneTwister, A::AbstractArray, g::RandIntGen)
+function rand!(mt::MersenneTwister, A::AbstractArray, g::RangeGenerator)
     for i = 1 : length(A)
         @inbounds A[i] = rand(mt, g)
     end
     return A
 end
 
-rand!{T<:Union(Signed,Unsigned,Bool,Char)}(mt::MersenneTwister, A::AbstractArray, r::UnitRange{T}) = rand!(mt, A, RandIntGen(r))
+rand!{T<:Union(Signed,Unsigned,BigInt,Bool,Char)}(mt::MersenneTwister, A::AbstractArray, r::UnitRange{T}) = rand!(mt, A, inrange(r))
 
 function rand!(mt::MersenneTwister, A::AbstractArray, r::AbstractArray)
-    g = RandIntGen(1:(length(r)))
+    g = inrange(1:(length(r)))
     for i = 1 : length(A)
         @inbounds A[i] = r[rand(mt, g)]
     end

test/random.jl

@@ -54,7 +54,7 @@ randn!(MersenneTwister(42), A)
 @test A == [-0.5560268761463861  0.027155338009193845;
             -0.444383357109696  -0.29948409035891055]
 
-for T in (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128,
+for T in (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, BigInt,
           Float16, Float32, Float64, Rational{Int})
     r = rand(convert(T, 97):convert(T, 122))
     @test typeof(r) == T
@@ -64,7 +64,7 @@ for T in (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt
     @test 97 <= r <= 122
     @test mod(r,2)==1
 
-    if T<:Integer
+    if T<:Integer && !(T===BigInt)
         x = rand(typemin(T):typemax(T))
         @test isa(x,T)
         @test typemin(T) <= x <= typemax(T)
@@ -75,11 +75,26 @@ 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.inrange(uint64(1:k)).u for k in 13 .+ int64(2).^(32:62)])
+    @test all([div(0x00010000000000000000,k)*k - 1 == Base.Random.inrange(int64(1:k)).u for k in 13 .+ int64(2).^(32:61)])
 
 end
 
+# BigInt specific
+for T in [UInt32, UInt64, UInt128, Int128]
+    s = big(typemax(T)-1000) : big(typemax(T)) + 10000
+    @test rand(s) != rand(s)
+    @test big(typemax(T)-1000) <= rand(s) <= big(typemax(T)) + 10000
+    r = rand(s, 1, 2)
+    @test size(r) == (1, 2)
+    @test typeof(r) == Matrix{BigInt}
+
+    srand(0)
+    r = rand(s)
+    srand(0)
+    @test rand(s) == r
+end
+
 randn(100000)
 randn!(Array(Float64, 100000))
 randn(MersenneTwister(10), 100000)
@@ -181,8 +196,8 @@ 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.inrange(uint64(1:k)).u for k in 13 .+ int64(2).^(1:30)])
+@test all([div(0x000100000000,k)*k - 1 == Base.Random.inrange(int64(1:k)).u for k in 13 .+ int64(2).^(1:30)])
 
 import Base.Random: uuid4, UUID