Implement rand()

Project.toml

@@ -5,6 +5,7 @@ version = "0.4.10"
 [deps]
 DecFP_jll = "47200ebd-12ce-5be5-abb7-8e082af23329"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]

src/DecFP.jl

@@ -3,7 +3,7 @@ module DecFP
 
 using DecFP_jll
 
-import Printf, SpecialFunctions
+import Printf, Random, SpecialFunctions
 
 export Dec32, Dec64, Dec128, @d_str, @d32_str, @d64_str, @d128_str, exponent10, ldexp10
 
@@ -497,6 +497,8 @@ for w in (32,64,128)
 
     @eval SpecialFunctions.gamma(x::$BID) = @xchk(ccall(($(bidsym(w,:tgamma)), libbid), $BID, ($BID,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)
 
+    @eval Random.rand(r::Random.AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{$BID}}) =  $BID(1, rand(r, zero($Ti):$Ti(maxintfloat($BID)) - one($Ti)), -(9 * $w ÷ 32 - 2))
+
     for (r,c) in ((RoundingMode{:Nearest},"round_integral_nearest_even"), (RoundingMode{:NearestTiesAway},"round_integral_nearest_away"), (RoundingMode{:ToZero},"round_integral_zero"), (RoundingMode{:Up},"round_integral_positive"), (RoundingMode{:Down},"round_integral_negative"))
         @eval Base.round(x::$BID, ::$r) = @xchk(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)
     end

test/runtests.jl

@@ -1,10 +1,18 @@
-using DecFP, Test, Printf, Base.MathConstants, SpecialFunctions
+using DecFP, Test, Printf, Random, Base.MathConstants, SpecialFunctions
 
 @test DecFP.flags[Threads.threadid()] == 0
 
 import DecFP.isnanstr
 @test isnanstr("nan") && isnanstr("  +NAN") && isnanstr("-NaN") && !isnanstr("nano")
 
+function hist(X, n)
+    v = zeros(Int, n)
+    for x in X
+        v[floor(Int, x*n) + 1] += 1
+    end
+    v
+end
+
 function testthreads(T, i, mode)
     @test @sprintf("%.0f", T(i)) == string(i)
     @test rounding(T) == RoundNearest
@@ -236,6 +244,27 @@ for T in (Dec32, Dec64, Dec128)
         @test gamma(T(x)) ≈ gamma(x)
     end
 
+    Random.seed!(1234)
+    @test rand(T) != rand(T)
+    @test 0 <= rand(T) < 1
+    r = rand(convert(T, 97):convert(T, 122))
+    @test typeof(r) == T
+    @test 97 <= r <= 122
+    r = rand(convert(T, 97):convert(T, 2):convert(T, 122), 2)[1]
+    @test typeof(r) == T
+    @test 97 <= r <= 122
+    @test mod(r, 2)==1
+    A = Vector{T}(undef, 16)
+    rand!(A)
+    @test all(0 .<= A .< 1)
+    # test uniform distribution
+    # array version
+    counts = hist(rand(T, 2000), 4)
+    @test minimum(counts) > 300 # should fail with proba < 1e-26
+    # scalar version
+    counts = hist([rand(T) for i in 1:2000], 4)
+    @test minimum(counts) > 300
+
     for c in (π, e, γ, catalan, φ)
         @test T(c) ≈ Float64(c)
     end