From 79e4104af7ba2b2d1ff991802978fd196ef55a05 Mon Sep 17 00:00:00 2001 From: Rafael Fourquet Date: Mon, 8 Sep 2014 21:02:34 +0530 Subject: [PATCH] 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 bf8c452112 are re-added. --- base/gmp.jl | 17 +++++++++- base/random.jl | 92 +++++++++++++++++++++++++++++++++++++++++++------- base/sysimg.jl | 6 ++-- test/random.jl | 27 +++++++++++---- 4 files changed, 121 insertions(+), 21 deletions(-) diff --git a/base/gmp.jl b/base/gmp.jl index cb0d50d2c64f6..ccc491657acf0 100644 --- a/base/gmp.jl +++ b/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) diff --git a/base/random.jl b/base/random.jl index 620ef46adef0d..025b4f7012bdc 100644 --- a/base/random.jl +++ 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)<> 32 == 0 x = rand(mt, UInt32) @@ -402,7 +437,7 @@ 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) @@ -410,23 +445,56 @@ 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)) +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 diff --git a/base/sysimg.jl b/base/sysimg.jl index 74d6ee5589042..67aa85ba7dd77 100644 --- a/base/sysimg.jl +++ b/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") diff --git a/test/random.jl b/test/random.jl index c42f9860f0b42..221118dcfded3 100644 --- a/test/random.jl +++ b/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