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. This is fixed here, but a mecanism handling
  not-implemented types should be added.
* RandIntGen is renamed to RangeGeneratorInt.
* A type RangeGeneratorBigInt similar to RangeGeneratorInt (both
  subtypes of RangeGenerator) is introduced to handle rand on
  BigInt ranges. RangeGenerator is the generic constructor of such
  objects, taking a range as parameter.
  Note: two different versions are implemented depending on the
  GMP version (it is faster with version 6 on big ranges).
* BigInt tests from commit bf8c452 are re-added.

Author: rfourquet

base/gmp.jl

@@ -18,10 +18,25 @@ else
 end
 typealias CdoubleMax Union(Float16, Float32, Float64)
 
+const GMP_VERSION = VersionNumber(bytestring(unsafe_load(cglobal((:__gmp_version, :libgmp), Ptr{Cchar}))))
+const GMP_BITS_PER_LIMB = Int(unsafe_load(cglobal((:__gmp_bits_per_limb, :libgmp), Cint)))
+
+# GMP's mp_limb_t is by default a typedef of `unsigned long`, but can also be configured to be either
+# `unsigned int` or `unsigned long long int`. The correct unsigned type is here named Limb, and must
+# be used whenever mp_limb_t is in the signature of ccall'ed GMP functions.
+if GMP_BITS_PER_LIMB == 32
+    typealias Limb UInt32
+elseif GMP_BITS_PER_LIMB == 64
+    typealias Limb UInt64
+else
+    error("GMP: cannot determine the type mp_limb_t (__gmp_bits_per_limb == $GMP_BITS_PER_LIMB)")
+end
+
+
 type BigInt <: Integer
     alloc::Cint
     size::Cint
-    d::Ptr{Culong}
+    d::Ptr{Limb}
     function BigInt()
         b = new(zero(Cint), zero(Cint), C_NULL)
         ccall((:__gmpz_init,:libgmp), Void, (Ptr{BigInt},), &b)

base/random.jl

@@ -1,6 +1,7 @@
 module Random
 
 using Base.dSFMT
+using Base.GMP: GMP_VERSION, Limb
 
 export srand,
        rand, rand!,
@@ -362,31 +363,65 @@ 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 RangeGeneratorInt{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
 end
 # generators with 32, 128 bits entropy
-RandIntGen{T, U<:Union(UInt32, UInt128)}(a::T, k::U) = RandIntGen{T, U}(a, k, maxmultiple(k))
+RangeGeneratorInt{T, U<:Union(UInt32, UInt128)}(a::T, k::U) = RangeGeneratorInt{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))
+RangeGeneratorInt{T}(a::T, k::UInt64) = RangeGeneratorInt{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))
+RangeGenerator{T<:Unsigned}(r::UnitRange{T}) = isempty(r) ? error("range must be non-empty") : RangeGeneratorInt(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 RangeGenerator(r::UnitRange{$T}) = isempty(r) ? error("range must be non-empty") : RangeGeneratorInt(first(r), convert($U, unsigned(last(r) - first(r)) + one($U))) # overflow ok
+end
+
+if GMP_VERSION.major >= 6
+    immutable RangeGeneratorBigInt <: RangeGenerator
+        a::BigInt             # first
+        m::BigInt             # range length - 1
+        nlimbs::Int           # number of limbs in generated BigInt's
+        mask::Limb            # applied to the highest limb
+    end
+
+else
+    immutable RangeGeneratorBigInt <: RangeGenerator
+        a::BigInt             # first
+        m::BigInt             # range length - 1
+        limbs::Array{Limb}    # buffer to be copied into generated BigInt's
+        mask::Limb            # applied to the highest limb
+
+        RangeGeneratorBigInt(a, m, nlimbs, mask) = new(a, m, Array(Limb, nlimbs), mask)
+    end
+end
+
+
+
+function RangeGenerator(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(Limb))
+    highbits > 0 && (nlimbs += 1)
+    mask = highbits == 0 ? ~zero(Limb) : one(Limb)<<highbits - one(Limb)
+    return RangeGeneratorBigInt(first(r), m, 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})
+# RangeGeneratorInt is responsible for providing the right value of k
+function rand{T<:Union(UInt64, Int64)}(mt::MersenneTwister, g::RangeGeneratorInt{T,UInt64})
     local x::UInt64
     if (g.k - 1) >> 32 == 0
         x = rand(mt, UInt32)
@@ -402,31 +437,64 @@ function rand{T<:Union(UInt64, Int64)}(mt::MersenneTwister, g::RandIntGen{T,UInt
     return reinterpret(T, reinterpret(UInt64, g.a) + rem_knuth(x, g.k))
 end
 
-function rand{T<:Integer, U<:Unsigned}(mt::MersenneTwister, g::RandIntGen{T,U})
+function rand{T<:Integer, U<:Unsigned}(mt::MersenneTwister, g::RangeGeneratorInt{T,U})
     x = rand(mt, U)
     while x > g.u
         x = rand(mt, U)
     end
     (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))
+if GMP_VERSION.major >= 6
+    # mpz_limbs_write and mpz_limbs_finish are available only in GMP version 6
+    function rand(rng::AbstractRNG, g::RangeGeneratorBigInt)
+        x = BigInt()
+        while true
+            # note: on CRAY computers, the second argument may be of type Cint (48 bits) and not Clong
+            xd = ccall((:__gmpz_limbs_write, :libgmp), Ptr{Limb}, (Ptr{BigInt}, Clong), &x, g.nlimbs)
+            limbs = pointer_to_array(xd, g.nlimbs)
+            rand!(rng, limbs)
+            limbs[end] &= g.mask
+            ccall((:__gmpz_limbs_finish, :libgmp), Void, (Ptr{BigInt}, Clong), &x, g.nlimbs)
+            x <= g.m && break
+        end
+        ccall((:__gmpz_add, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &x, &x, &g.a)
+        return x
+    end
+else
+    function rand(rng::AbstractRNG, g::RangeGeneratorBigInt)
+        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(Limb), 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
+end
+
+rand{T<:Union(Signed,Unsigned,BigInt,Bool,Char)}(mt::MersenneTwister, r::UnitRange{T}) = rand(mt, RangeGenerator(r))
+
 
 # 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, RangeGenerator(r))
 
 function rand!(mt::MersenneTwister, A::AbstractArray, r::AbstractArray)
-    g = RandIntGen(1:(length(r)))
+    g = RangeGenerator(1:(length(r)))
     for i = 1 : length(A)
         @inbounds A[i] = r[rand(mt, g)]
     end

base/sysimg.jl

@@ -176,6 +176,9 @@ include("combinatorics.jl")
 include("rounding.jl")
 importall .Rounding
 
+# version
+include("version.jl")
+
 # BigInts and BigFloats
 include("gmp.jl")
 importall .GMP
@@ -198,8 +201,7 @@ include("darray.jl")
 include("mmap.jl")
 include("sharedarray.jl")
 
-# utilities - version, timing, help, edit, metaprogramming
-include("version.jl")
+# utilities - timing, help, edit, metaprogramming
 include("datafmt.jl")
 importall .DataFmt
 include("deepcopy.jl")

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.RangeGenerator(uint64(1:k)).u for k in 13 .+ int64(2).^(32:62)])
+    @test all([div(0x00010000000000000000,k)*k - 1 == Base.Random.RangeGenerator(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()
 randn(100000)
 randn!(Array(Float64, 100000))
@@ -190,8 +205,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.RangeGenerator(uint64(1:k)).u for k in 13 .+ int64(2).^(1:30)])
+@test all([div(0x000100000000,k)*k - 1 == Base.Random.RangeGenerator(int64(1:k)).u for k in 13 .+ int64(2).^(1:30)])
 
 import Base.Random: uuid4, UUID