Merge pull request #22720 from JuliaLang/rf/rand-BigFloat

implement rand(::BigFloat) (close #13948)

NEWS.md

@@ -208,6 +208,8 @@ Library improvements
   * Mutating versions of `randperm` and `randcycle` have been added:
     `randperm!` and `randcycle!` ([#22723]).
 
+  * `BigFloat` random numbers can now be generated ([#22720]).
+
 Compiler/Runtime improvements
 -----------------------------
 

base/random/RNGs.jl

@@ -45,12 +45,7 @@ RandomDevice
 
 ### generation of floats
 
-rand(rng::RandomDevice, ::Type{Close1Open2}) =
-    reinterpret(Float64, 0x3ff0000000000000 | rand(rng, UInt64) & 0x000fffffffffffff)
-
-rand(rng::RandomDevice, ::Type{CloseOpen}) = rand(rng, Close1Open2) - 1.0
-
-@inline rand(r::RandomDevice, T::BitFloatType) = rand_generic(r, T)
+@inline rand(r::RandomDevice, I::FloatInterval) = rand_generic(r, I)
 
 
 ## MersenneTwister
@@ -198,13 +193,13 @@ const GLOBAL_RNG = MersenneTwister(0)
 #### helper functions
 
 # precondition: !mt_empty(r)
-@inline rand_inbounds(r::MersenneTwister, ::Type{Close1Open2}) = mt_pop!(r)
-@inline rand_inbounds(r::MersenneTwister, ::Type{CloseOpen}) =
-    rand_inbounds(r, Close1Open2) - 1.0
-@inline rand_inbounds(r::MersenneTwister) = rand_inbounds(r, CloseOpen)
+@inline rand_inbounds(r::MersenneTwister, ::Close1Open2_64) = mt_pop!(r)
+@inline rand_inbounds(r::MersenneTwister, ::CloseOpen_64) =
+    rand_inbounds(r, Close1Open2()) - 1.0
+@inline rand_inbounds(r::MersenneTwister) = rand_inbounds(r, CloseOpen())
 
 @inline rand_ui52_raw_inbounds(r::MersenneTwister) =
-    reinterpret(UInt64, rand_inbounds(r, Close1Open2))
+    reinterpret(UInt64, rand_inbounds(r, Close1Open2()))
 @inline rand_ui52_raw(r::MersenneTwister) = (reserve_1(r); rand_ui52_raw_inbounds(r))
 
 @inline function rand_ui2x52_raw(r::MersenneTwister)
@@ -222,10 +217,9 @@ rand_ui23_raw(r::MersenneTwister) = rand_ui52_raw(r)
 
 #### floats
 
-@inline rand(r::MersenneTwister, ::Type{I}) where {I<:FloatInterval} =
-    (reserve_1(r); rand_inbounds(r, I))
+@inline rand(r::MersenneTwister, I::FloatInterval_64) = (reserve_1(r); rand_inbounds(r, I))
 
-@inline rand(r::MersenneTwister, T::BitFloatType) = rand_generic(r, T)
+@inline rand(r::MersenneTwister, I::FloatInterval) = rand_generic(r, I)
 
 #### integers
 
@@ -251,8 +245,7 @@ rand(r::MersenneTwister, ::Type{Int128}) = reinterpret(Int128, rand(r, UInt128))
 #### arrays of floats
 
 function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float64},
-                                     n=length(A),
-                                     ::Type{I}=CloseOpen) where I<:FloatInterval
+                                     n=length(A), I::FloatInterval_64=CloseOpen())
     # what follows is equivalent to this simple loop but more efficient:
     # for i=1:n
     #     @inbounds A[i] = rand(r, I)
@@ -275,14 +268,14 @@ end
 
 rand!(r::MersenneTwister, A::AbstractArray{Float64}) = rand_AbstractArray_Float64!(r, A)
 
-fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::Type{CloseOpen}) =
+fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::CloseOpen_64) =
     dsfmt_fill_array_close_open!(s, A, n)
 
-fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::Type{Close1Open2}) =
+fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::Close1Open2_64) =
     dsfmt_fill_array_close1_open2!(s, A, n)
 
 function rand!(r::MersenneTwister, A::Array{Float64}, n::Int=length(A),
-               ::Type{I}=CloseOpen) where I<:FloatInterval
+               I::FloatInterval_64=CloseOpen())
     # depending on the alignment of A, the data written by fill_array! may have
     # to be left-shifted by up to 15 bytes (cf. unsafe_copy! below) for
     # reproducibility purposes;
@@ -319,12 +312,12 @@ end
     (u & 0x007fffff007fffff007fffff007fffff) | 0x3f8000003f8000003f8000003f800000
 
 function rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}},
-               ::Type{Close1Open2})
+               ::Close1Open2_64)
     T = eltype(A)
     n = length(A)
     n128 = n * sizeof(T) ÷ 16
     rand!(r, unsafe_wrap(Array, convert(Ptr{Float64}, pointer(A)), 2*n128),
-          2*n128, Close1Open2)
+          2*n128, Close1Open2())
     A128 = unsafe_wrap(Array, convert(Ptr{UInt128}, pointer(A)), n128)
     @inbounds for i in 1:n128
         u = A128[i]
@@ -347,17 +340,17 @@ function rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}},
     A
 end
 
-function rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}},
-               ::Type{CloseOpen})
-    rand!(r, A, Close1Open2)
+function rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}}, ::CloseOpen_64)
+    rand!(r, A, Close1Open2())
     I32 = one(Float32)
     for i in eachindex(A)
         @inbounds A[i] = Float32(A[i])-I32 # faster than "A[i] -= one(T)" for T==Float16
     end
     A
 end
 
-rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}}) = rand!(r, A, CloseOpen)
+rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}}) =
+    rand!(r, A, CloseOpen())
 
 #### arrays of integers
 
@@ -368,7 +361,7 @@ function rand!(r::MersenneTwister, A::Array{UInt128}, n::Int=length(A))
     Af = unsafe_wrap(Array, convert(Ptr{Float64}, pointer(A)), 2n)
     i = n
     while true
-        rand!(r, Af, 2i, Close1Open2)
+        rand!(r, Af, 2i, Close1Open2())
         n < 5 && break
         i = 0
         @inbounds while n-i >= 5

base/random/generation.jl

@@ -10,25 +10,92 @@
 
 ### random floats
 
-@inline rand(r::AbstractRNG=GLOBAL_RNG) = rand(r, CloseOpen)
+# CloseOpen(T) is the fallback for an AbstractFloat T
+@inline rand(r::AbstractRNG=GLOBAL_RNG, ::Type{T}=Float64) where {T<:AbstractFloat} =
+    rand(r, CloseOpen(T))
 
 # generic random generation function which can be used by RNG implementors
 # it is not defined as a fallback rand method as this could create ambiguities
-@inline rand_generic(r::AbstractRNG, ::Type{Float64}) = rand(r, CloseOpen)
 
-rand_generic(r::AbstractRNG, ::Type{Float16}) =
+rand_generic(r::AbstractRNG, ::CloseOpen{Float16}) =
     Float16(reinterpret(Float32,
                         (rand_ui10_raw(r) % UInt32 << 13) & 0x007fe000 | 0x3f800000) - 1)
 
-rand_generic(r::AbstractRNG, ::Type{Float32}) =
+rand_generic(r::AbstractRNG, ::CloseOpen{Float32}) =
     reinterpret(Float32, rand_ui23_raw(r) % UInt32 & 0x007fffff | 0x3f800000) - 1
 
+rand_generic(r::AbstractRNG, ::Close1Open2_64) =
+    reinterpret(Float64, 0x3ff0000000000000 | rand(r, UInt64) & 0x000fffffffffffff)
+
+rand_generic(r::AbstractRNG, ::CloseOpen_64) = rand(r, Close1Open2()) - 1.0
+
+#### BigFloat
+
+const bits_in_Limb = sizeof(Limb) << 3
+const Limb_high_bit = one(Limb) << (bits_in_Limb-1)
+
+struct BigFloatRandGenerator
+    prec::Int
+    nlimbs::Int
+    limbs::Vector{Limb}
+    shift::UInt
+
+    function BigFloatRandGenerator(prec::Int=precision(BigFloat))
+        nlimbs = (prec-1) ÷ bits_in_Limb + 1
+        limbs = Vector{Limb}(nlimbs)
+        shift = nlimbs * bits_in_Limb - prec
+        new(prec, nlimbs, limbs, shift)
+    end
+end
+
+function _rand(rng::AbstractRNG, gen::BigFloatRandGenerator)
+    z = BigFloat()
+    limbs = gen.limbs
+    rand!(rng, limbs)
+    @inbounds begin
+        limbs[1] <<= gen.shift
+        randbool = iszero(limbs[end] & Limb_high_bit)
+        limbs[end] |= Limb_high_bit
+    end
+    z.sign = 1
+    unsafe_copy!(z.d, pointer(limbs), gen.nlimbs)
+    (z, randbool)
+end
+
+function rand(rng::AbstractRNG, gen::BigFloatRandGenerator, ::Close1Open2{BigFloat})
+    z = _rand(rng, gen)[1]
+    z.exp = 1
+    z
+end
+
+function rand(rng::AbstractRNG, gen::BigFloatRandGenerator, ::CloseOpen{BigFloat})
+    z, randbool = _rand(rng, gen)
+    z.exp = 0
+    randbool &&
+        ccall((:mpfr_sub_d, :libmpfr), Int32,
+              (Ptr{BigFloat}, Ptr{BigFloat}, Cdouble, Int32),
+              &z, &z, 0.5, Base.MPFR.ROUNDING_MODE[])
+    z
+end
+
+# alternative, with 1 bit less of precision
+# TODO: make an API for requesting full or not-full precision
+function rand(rng::AbstractRNG, gen::BigFloatRandGenerator, ::CloseOpen{BigFloat}, ::Void)
+    z = rand(rng, Close1Open2(BigFloat), gen)
+    ccall((:mpfr_sub_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32),
+          &z, &z, 1, Base.MPFR.ROUNDING_MODE[])
+    z
+end
+
+rand_generic(rng::AbstractRNG, I::FloatInterval{BigFloat}) =
+    rand(rng, BigFloatRandGenerator(), I)
+
 ### random integers
 
 rand_ui10_raw(r::AbstractRNG) = rand(r, UInt16)
 rand_ui23_raw(r::AbstractRNG) = rand(r, UInt32)
 
-@inline rand_ui52_raw(r::AbstractRNG) = reinterpret(UInt64, rand(r, Close1Open2))
+@inline rand_ui52_raw(r::AbstractRNG) = reinterpret(UInt64, rand(r, Close1Open2()))
 @inline rand_ui52(r::AbstractRNG) = rand_ui52_raw(r) & 0x000fffffffffffff
 
 ### random complex numbers
@@ -72,6 +139,27 @@ rand(                T::Type, d::Integer, dims::Integer...) = rand(T, Dims((d, d
 # rand(r, ()) would match both this method and rand(r, dims::Dims)
 # moreover, a call like rand(r, NotImplementedType()) would be an infinite loop
 
+#### arrays of floats
+
+rand!(r::AbstractRNG, A::AbstractArray, ::Type{T}) where {T<:AbstractFloat} =
+    rand!(r, A, CloseOpen{T}())
+
+function rand!(r::AbstractRNG, A::AbstractArray, I::FloatInterval)
+    for i in eachindex(A)
+        @inbounds A[i] = rand(r, I)
+    end
+    A
+end
+
+function rand!(rng::AbstractRNG, A::AbstractArray, I::FloatInterval{BigFloat})
+    gen = BigFloatRandGenerator()
+    for i in eachindex(A)
+        @inbounds A[i] = rand(rng, gen, I)
+    end
+    A
+end
+
+rand!(A::AbstractArray, I::FloatInterval) = rand!(GLOBAL_RNG, A, I)
 
 ## Generate random integer within a range
 

base/random/random.jl

@@ -22,9 +22,17 @@ export srand,
 
 abstract type AbstractRNG end
 
-abstract type FloatInterval end
-mutable struct CloseOpen <: FloatInterval end
-mutable struct Close1Open2 <: FloatInterval end
+abstract type FloatInterval{T<:AbstractFloat} end
+
+struct CloseOpen{  T<:AbstractFloat} <: FloatInterval{T} end # interval [0,1)
+struct Close1Open2{T<:AbstractFloat} <: FloatInterval{T} end # interval [1,2)
+
+const FloatInterval_64 = FloatInterval{Float64}
+const CloseOpen_64     = CloseOpen{Float64}
+const Close1Open2_64   = Close1Open2{Float64}
+
+CloseOpen(  ::Type{T}=Float64) where {T<:AbstractFloat} = CloseOpen{T}()
+Close1Open2(::Type{T}=Float64) where {T<:AbstractFloat} = Close1Open2{T}()
 
 const BitFloatType = Union{Type{Float16},Type{Float32},Type{Float64}}
 

doc/src/stdlib/numbers.md

@@ -144,12 +144,12 @@ dimension specifications `dims...` (which can be given as a tuple) to generate a
 values.
 
 A `MersenneTwister` or `RandomDevice` RNG can generate random numbers of the following types:
-[`Float16`](@ref), [`Float32`](@ref), [`Float64`](@ref), [`Bool`](@ref), [`Int8`](@ref),
-[`UInt8`](@ref), [`Int16`](@ref), [`UInt16`](@ref), [`Int32`](@ref), [`UInt32`](@ref),
-[`Int64`](@ref), [`UInt64`](@ref), [`Int128`](@ref), [`UInt128`](@ref), [`BigInt`](@ref)
-(or complex numbers of those types). Random floating point numbers are generated uniformly
-in ``[0, 1)``. As `BigInt` represents unbounded integers, the interval must be specified
-(e.g. `rand(big(1:6))`).
+[`Float16`](@ref), [`Float32`](@ref), [`Float64`](@ref), [`BigFloat`](@ref), [`Bool`](@ref),
+[`Int8`](@ref), [`UInt8`](@ref), [`Int16`](@ref), [`UInt16`](@ref), [`Int32`](@ref),
+[`UInt32`](@ref), [`Int64`](@ref), [`UInt64`](@ref), [`Int128`](@ref), [`UInt128`](@ref),
+[`BigInt`](@ref) (or complex numbers of those types).
+Random floating point numbers are generated uniformly in ``[0, 1)``. As `BigInt` represents
+unbounded integers, the interval must be specified (e.g. `rand(big(1:6))`).
 
 ```@docs
 Base.Random.srand

test/random.jl

@@ -308,7 +308,7 @@ end
 for rng in ([], [MersenneTwister(0)], [RandomDevice()])
     ftypes = [Float16, Float32, Float64]
     cftypes = [Complex32, Complex64, Complex128, ftypes...]
-    types = [Bool, Char, Base.BitInteger_types..., ftypes...]
+    types = [Bool, Char, BigFloat, Base.BitInteger_types..., ftypes...]
     randset = Set(rand(Int, 20))
     randdict = Dict(zip(rand(Int,10), rand(Int, 10)))
     collections = [IntSet(rand(1:100, 20)) => Int,
@@ -322,6 +322,9 @@ for rng in ([], [MersenneTwister(0)], [RandomDevice()])
                    Float64                 => Float64,
                    "qwèrtï"                => Char,
                    GenericString("qwèrtï") => Char]
+    functypes = Dict(rand  => types, randn  => cftypes, randexp  => ftypes,
+                     rand! => types, randn! => cftypes, randexp! => ftypes)
+
     b2 = big(2)
     u3 = UInt(3)
     for f in [rand, randn, randexp]
@@ -330,7 +333,7 @@ for rng in ([], [MersenneTwister(0)], [RandomDevice()])
         f(rng..., 2, 3)               ::Array{Float64, 2}
         f(rng..., b2, u3)             ::Array{Float64, 2}
         f(rng..., (2, 3))             ::Array{Float64, 2}
-        for T in (f === rand ? types : f === randn ? cftypes : ftypes)
+        for T in functypes[f]
             a0 = f(rng..., T)         ::T
             a1 = f(rng..., T, 5)      ::Vector{T}
             a2 = f(rng..., T, 2, 3)   ::Array{T, 2}
@@ -370,7 +373,7 @@ for rng in ([], [MersenneTwister(0)], [RandomDevice()])
         @test_throws ArgumentError rand(rng..., C, 5)
     end
     for f! in [rand!, randn!, randexp!]
-        for T in (f! === rand! ? types : f! === randn! ? cftypes : ftypes)
+        for T in functypes[f!]
             X = T == Bool ? T[0,1] : T[0,1,2]
             for A in (Array{T}(5), Array{T}(2, 3), GenericArray{T}(5), GenericArray{T}(2, 3))
                 f!(rng..., A)                    ::typeof(A)
@@ -409,17 +412,19 @@ for rng in ([], [MersenneTwister(0)], [RandomDevice()])
     end
 end
 
-function hist(X,n)
-    v = zeros(Int,n)
+function hist(X, n)
+    v = zeros(Int, n)
     for x in X
-        v[floor(Int,x*n)+1] += 1
+        v[floor(Int, x*n) + 1] += 1
     end
     v
 end
 
 # test uniform distribution of floats
-for rng in [srand(MersenneTwister(0)), RandomDevice()]
-    for T in [Float16,Float32,Float64]
+for rng in [srand(MersenneTwister(0)), RandomDevice()],
+    T in [Float16, Float32, Float64, BigFloat],
+        prec in (T == BigFloat ? [3, 53, 64, 100, 256, 1000] : [256])
+    setprecision(BigFloat, prec) do
         # array version
         counts = hist(rand(rng, T, 2000), 4)
         @test minimum(counts) > 300 # should fail with proba < 1e-26