rand BFloat16 sampling without conversion

src/bfloat16.jl

@@ -350,7 +350,17 @@ Printf.tofloat(x::BFloat16) = Float32(x)
 
 # Random
 import Random: rand, randn, randexp, AbstractRNG, Sampler
-rand(rng::AbstractRNG, ::Sampler{BFloat16}) = convert(BFloat16, rand(rng))
+
+"""Sample a BFloat16 from [0,1) by setting random mantissa
+bits for one(BFloat16) to obtain [1,2) (where floats are uniformly
+distributed) then subtract 1 for [0,1)."""
+function rand(rng::AbstractRNG, ::Sampler{BFloat16})
+    u = reinterpret(UInt16, one(BFloat16))
+    # shift random bits into BFloat16 mantissa (1 sign + 8 exp bits = 9)
+    u |= rand(rng, UInt16) >> 9                     # u in [1,2)
+    return reinterpret(BFloat16, u) - one(BFloat16) # -1 for [0,1)
+end
+
 randn(rng::AbstractRNG, ::Type{BFloat16}) = convert(BFloat16, randn(rng))
 randexp(rng::AbstractRNG, ::Type{BFloat16}) = convert(BFloat16, randexp(rng))
 

test/runtests.jl

@@ -175,5 +175,18 @@ end
   @test a-1+1 == a      # but -1 can
 end
 
+@testset "rand sampling" begin
+  Random.seed(123)
+  mi, ma = extrema(rand(BFloat16, 1_000_000))
+
+  # zero should be the lowest BFloat16 sampled
+  @test mi === zero(BFloat16)
+
+  # prevfloat(one(BFloat16)) cannot be sampled bc
+  # prevfloat(BFloat16(2)) - 1 is _two_ before one(BFloat16)
+  # (a statistical flaw of the [1,2)-1 sampling)
+  @test ma === prevfloat(one(BFloat16), 2)
+end
+
 include("structure.jl")
-include("mathfuncs.jl")
+include("mathfuncs.jl")