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

Author: rfourquet

NEWS.md

@@ -200,6 +200,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/generation.jl

@@ -29,6 +29,67 @@ rand_generic(r::AbstractRNG, ::Close1Open2_64) =
 
 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_generic(rng::AbstractRNG, ::Close1Open2{BigFloat},
+                      gen::BigFloatRandGenerator=BigFloatRandGenerator())
+    z = _rand(rng, gen)[1]
+    z.exp = 1
+    z
+end
+
+function rand_generic(rng::AbstractRNG, ::CloseOpen{BigFloat},
+                      gen::BigFloatRandGenerator=BigFloatRandGenerator())
+    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_generic(rng::AbstractRNG, ::CloseOpen{BigFloat}, ::Val{false},
+                      gen::BigFloatRandGenerator=BigFloatRandGenerator())
+    z = rand_generic(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
+
 ### random integers
 
 rand_ui10_raw(r::AbstractRNG) = rand(r, UInt16)
@@ -90,6 +151,14 @@ function rand!(r::AbstractRNG, A::AbstractArray, I::FloatInterval)
     A
 end
 
+function rand!(rng::AbstractRNG, A::AbstractArray, I::FloatInterval{BigFloat})
+    gen = BigFloatRandGenerator()
+    for i in eachindex(A)
+        @inbounds A[i] = rand_generic(rng, I, gen)
+    end
+    A
+end
+
 rand!(A::AbstractArray, I::FloatInterval) = rand!(GLOBAL_RNG, A, I)
 
 ## Generate random integer within a range

doc/src/stdlib/numbers.md

@@ -143,12 +143,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), [`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))`).
 
 ```@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