Skip to content

Commit

Permalink
Merge pull request #9122 from JuliaLang/rf/rand-bigrange
Browse files Browse the repository at this point in the history
Implement rand(::Range{BigInt})
  • Loading branch information
ViralBShah committed Dec 2, 2014
2 parents bf57363 + 79e4104 commit 7eb9a14
Show file tree
Hide file tree
Showing 4 changed files with 121 additions and 21 deletions.
17 changes: 16 additions & 1 deletion base/gmp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
92 changes: 80 additions & 12 deletions base/random.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
module Random

using Base.dSFMT
using Base.GMP: GMP_VERSION, Limb

export srand,
rand, rand!,
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
6 changes: 4 additions & 2 deletions base/sysimg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,9 @@ include("combinatorics.jl")
include("rounding.jl")
importall .Rounding

# version
include("version.jl")

# BigInts and BigFloats
include("gmp.jl")
importall .GMP
Expand All @@ -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")
Expand Down
27 changes: 21 additions & 6 deletions test/random.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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))
Expand Down Expand Up @@ -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

Expand Down

0 comments on commit 7eb9a14

Please sign in to comment.