Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor rand on UnitRange / AbstractArray #9324

Closed
wants to merge 1 commit into from
Closed

refactor rand on UnitRange / AbstractArray #9324

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 51 additions & 73 deletions base/random.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
module Random

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

export srand,
Expand Down Expand Up @@ -403,97 +403,85 @@ rem_knuth{T<:Unsigned}(a::T, b::T) = b != 0 ? a % b : a
# 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) = (div(0x0000000100000000,k + (k == 0))*k - 1) % UInt32
maxmultiple(k::UInt64) = (div(0x00000000000000010000000000000000, k + (k == 0))*k - 1) % UInt64

maxmultiple64(k::UInt64) = (div(0x00000000000000010000000000000000, k + (k == 0))*k - 1) % UInt64

# maximum multiple of k within 1:2^32 or 1:2^64, depending on size
maxmultiple(k::UInt64) = (div((k >> 32 != 0)*0x0000000000000000FFFFFFFF00000000 + 0x0000000100000000, k + (k == 0))*k - 1) % UInt64

# 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) = (div((k >> 32 != 0)*0x0000000000000000FFFFFFFF00000000 + 0x0000000100000000, k + (k == 0))*k - 1) % UInt64

abstract RangeGenerator
# utility to generate random numbers in a range
abstract RangeGenerator{T<:Integer}

immutable RangeGeneratorInt{T<:Integer, U<:Unsigned} <: RangeGenerator
a::T # first element of the range
immutable RangeGeneratorInt{T<:Integer, U<:Union(UInt32, UInt64, UInt128)} <: RangeGenerator{T}
k::U # range length or zero for full range
u::U # rejection threshold
end
# generators with 32, 128 bits entropy
RangeGeneratorInt{T, U<:Union(UInt32, UInt128)}(a::T, k::U) = RangeGeneratorInt{T, U}(a, k, maxmultiple(k))
# mixed 32/64 bits entropy generator
RangeGeneratorInt{T}(a::T, k::UInt64) = RangeGeneratorInt{T,UInt64}(a, k, maxmultiplemix(k))


# generator for ranges
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 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
@inline RangeGenerator{T, U}(::Type{T}, m::U) = (k=m+one(U); RangeGeneratorInt{T, U}(k, maxmultiple(k)))

if GMP_VERSION.major >= 6
immutable RangeGeneratorBigInt <: RangeGenerator
a::BigInt # first
immutable RangeGeneratorBigInt <: RangeGenerator{BigInt}
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
immutable RangeGeneratorBigInt <: RangeGenerator{BigInt}
m::BigInt # range length - 1
limbs::Vector{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)
RangeGeneratorBigInt(m, nlimbs, mask) = new(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")
function RangeGenerator(::Type{BigInt}, m::BigInt)
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)
return RangeGeneratorBigInt(m, nlimbs, mask)
end

# RangeGenerator dispatch

# this function uses 32 bit entropy for small ranges of length <= typemax(UInt32) + 1
# RangeGeneratorInt is responsible for providing the right value of k
function rand{T<:Union(UInt64, Int64)}(rng::AbstractRNG, g::RangeGeneratorInt{T,UInt64})
local x::UInt64
if (g.k - 1) >> 32 == 0
x = rand(rng, UInt32)
while x > g.u
x = rand(rng, UInt32)
end
else
x = rand(rng, UInt64)
while x > g.u
x = rand(rng, UInt64)
end
end
return reinterpret(T, reinterpret(UInt64, g.a) + rem_knuth(x, g.k))
end
@inline checknonempty(r) = isempty(r) && error("collection must be non-empty")
@inline RangeGenerator(r::AbstractArray) = (checknonempty(r); n = length(r); RangeGenerator(n-one(n)))
@inline RangeGenerator{T<:Union(IntTypes...,BigInt)}(r::UnitRange{T}) = (checknonempty(r); RangeGenerator(last(r)-first(r)))

@inline RangeGenerator{T<:Union(UInt32,UInt64,UInt128,BigInt)}(m::T) = RangeGenerator(T, m)
@inline RangeGenerator{T<:Union(Int32,Int64,Int128)}(m::T) = RangeGenerator(T, unsigned(m))
@inline RangeGenerator{T<:Union(Int8,UInt8,Int16,UInt16)}(m::T) = RangeGenerator(T, m % UInt32)

# rand on RangeGenerator

function rand{T<:Integer, U<:Unsigned}(rng::AbstractRNG, g::RangeGeneratorInt{T,U})
x = rand(rng, U)
while x > g.u
@inline rand{T}(rng::AbstractRNG, g::RangeGenerator{T}) = rand(rng, g, one(T))

@inline function rand_lessthan{U}(rng::AbstractRNG, u::U)
while true
x = rand(rng, U)
x <= u && return x
end
(unsigned(g.a) + rem_knuth(x, g.k)) % T
end

# this function uses 32 bit entropy for small ranges of length <= typemax(UInt32) + 1
# RangeGeneratorInt is responsible for providing the right value of k
@inline rand{T}(rng::AbstractRNG, g::RangeGeneratorInt{T,UInt64}, start::T) =
(g.k - 1) >> 32 == 0 ?
start + rand_lessthan(rng, g.u % UInt32) % UInt64 % g.k % T :
rand(rng, g, start, true)

@inline rand{T}(rng::AbstractRNG, g::RangeGeneratorInt{T}, start::T, dummy=true) =
start + rem_knuth(rand_lessthan(rng, g.u), g.k) % T


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)
function rand(rng::AbstractRNG, g::RangeGeneratorBigInt, start::BigInt)
x = BigInt()
while true
# note: on CRAY computers, the second argument may be of type Cint (48 bits) and not Clong
Expand All @@ -504,11 +492,11 @@ if GMP_VERSION.major >= 6
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)
ccall((:__gmpz_add, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &x, &x, &start)
return x
end
else
function rand(rng::AbstractRNG, g::RangeGeneratorBigInt)
function rand(rng::AbstractRNG, g::RangeGeneratorBigInt, start::BigInt)
x = BigInt()
while true
rand!(rng, g.limbs)
Expand All @@ -518,31 +506,21 @@ else
&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)
ccall((:__gmpz_add, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &x, &x, &start)
return x
end
end

rand{T<:Union(Signed,Unsigned,BigInt,Bool,Char)}(rng::AbstractRNG, r::UnitRange{T}) = rand(rng, 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(rng::AbstractRNG, r::AbstractArray) = @inbounds return r[rand(rng, 1:length(r))]

function rand!(rng::AbstractRNG, A::AbstractArray, g::RangeGenerator)
for i = 1 : length(A)
@inbounds A[i] = rand(rng, g)
end
return A
end
@inline rand(rng::AbstractRNG, r::AbstractArray) = rand(rng, r, RangeGenerator(r))

rand!{T<:Union(Signed,Unsigned,BigInt,Bool,Char)}(rng::AbstractRNG, A::AbstractArray, r::UnitRange{T}) = rand!(rng, A, RangeGenerator(r))
@inline rand(rng::AbstractRNG, r::AbstractArray, g::RangeGenerator) = @inbounds return r[rand(rng, g)]
@inline rand{T<:Union(IntTypes...,BigInt)}(rng::AbstractRNG, r::UnitRange{T}, g::RangeGenerator) = rand(rng, g, first(r))

function rand!(rng::AbstractRNG, A::AbstractArray, r::AbstractArray)
g = RangeGenerator(1:(length(r)))
function rand!(rng::AbstractRNG, A::AbstractArray, r::AbstractArray, g::RangeGenerator=RangeGenerator(r))
for i = 1 : length(A)
@inbounds A[i] = r[rand(rng, g)]
@inbounds A[i] = rand(rng, r, g)
end
return A
end
Expand Down