# Output distribution of rand(Float32) and rand(Float64), thread 2

**URL:** https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184
**Category:** Internals & Design
**Tags:** float, random
**Created:** [October 19, 2023, 1:10pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184 "2023-10-19T13:10:07Z")
**Posts on this page:** 1
**Showing post:** 261

<div class="post-metadata">

### Author: ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)
#### Post date: [October 30, 2023, 5:36pm UTC](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184/261 "2023-10-30T17:36:03Z")

</div>

Here’s the magic to jam @mikmoore’s latest Float32 N=32 version into XoshiroSimd, using a very naive copy-the-code\_llvm approach:

> **magic**
>
> Here’s the complete repl session:
> 
> ```Julia-repl
> ▶ ./julia
> _
> _ _ _(_)_ | Documentation: https://docs.julialang.org
> (_) | (_) (_) |
> _ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
> | | | | | | |/ _` | |
> | | |_| | | | (_| | | Version 1.11.0-DEV.595 (2023-10-03)
> _/ |\ __'_|_|_|\__'_| | Commit b790cf8f0c* (26 days old master)
> |__/ |
> 
> julia> using Random, BenchmarkTools
> 
> julia> A = Array{Float32}(undef, 2^24);
> 
> julia> @benchmark rand!(A)
> BenchmarkTools.Trial: 949 samples with 1 evaluation.
> Range (min … max): 5.200 ms … 6.274 ms ┊ GC (min … max): 0.00% … 0.00%
> Time (median): 5.240 ms ┊ GC (median): 0.00%
> Time (mean ± σ): 5.268 ms ± 83.680 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
> 
> ▃█▅▃▁▂
> ▄▅▅███████▆▇▆▆▆▇▅▆▆▆▄▆▅▇▆▆▅▄▄▅▁▅▄▄▅▄▁▄▄▁▁▄▄▄▁▄▁▁▁▁▄▁▄▁▁▄▁▄ ▇
> 5.2 ms Histogram: log(frequency) by time 5.7 ms <
> 
> Memory estimate: 0 bytes, allocs estimate: 0.
> 
> julia> @eval Random.XoshiroSimd begin
> @inline _bits2float(x::UInt64, ::Type{Float32}) = _f32bits(x)
> 
> @inline function _f32bits(u::UInt32)
> # map UIntXX linearly into FloatXX on range [0,1)
> U = UInt32
> F = Base.floattype(U)
> bitwidth = 8 * sizeof(U) # bits in type U
> lomask = typemax(U) >>> (bitwidth ÷ 2) # mask of low bits of u
> scale = F(exp2(-bitwidth))
> lo = F(u & lomask) * scale # must scale here for UInt16, because Float16(typemax(UInt16)) overflows
> hi = F(u & ~lomask) * scale # must scale here for UInt16, because Float16(typemax(UInt16)) overflows
> f0 = lo + hi # full-width value
> efflo = f0 - hi # get effective value of lo used
> roundedup = efflo > lo # determine whether f0 was rounded upward
> fu = reinterpret(U, f0)
> fu -= roundedup # if we rounded up, decrement to round down instead
> return fu
> end
> @inline function _f32bits(x::UInt64)
> ui = (x>>>32) % UInt32
> li = x % UInt32
> return (UInt64(_f32bits(ui)) << 32) | UInt64(_f32bits(li))
> end
> @inline _f32bits(x::VecElement{UInt64}) = _f32bits(x.value)
> end
> for N in [4,8,16]
> VT = :(NTuple{$N, VecElement{UInt64}})
> code = """
> %a = bitcast <$N x i64> %0 to <$(2N) x i32>
> %lowmask = shufflevector <1 x i32> <i32 65535>, <1 x i32> undef, <$(2N) x i32> zeroinitializer
> %b = and <$(2N) x i32> %a, %lowmask
> %c = uitofp <$(2N) x i32> %b to <$(2N) x float>
> %scale = shufflevector <1 x float> <float 0x3DF0000000000000>, <1 x float> undef, <$(2N) x i32> zeroinitializer
> %d = fmul <$(2N) x float> %c, %scale
> %himask = shufflevector <1 x i32> <i32 -65536>, <1 x i32> undef, <$(2N) x i32> zeroinitializer
> %e = and <$(2N) x i32> %a, %himask
> %f = uitofp <$(2N) x i32> %e to <$(2N) x float>
> %g = fmul <$(2N) x float> %f, %scale
> %h = fadd <$(2N) x float> %d, %g
> %i = fsub <$(2N) x float> %h, %g
> %j = fcmp olt <$(2N) x float> %d, %i
> %k = bitcast <$(2N) x float> %h to <$(2N) x i32>
> %l = sext <$(2N) x i1> %j to <$(2N) x i32>
> %m = add <$(2N) x i32> %l, %k
> %n = bitcast <$(2N) x i32> %m to <$N x i64>
> ret <$N x i64> %n
> """
> @eval Random.XoshiroSimd @inline _bits2float(x::$VT, ::Type{Float32}) = llvmcall($code, $VT, Tuple{$VT}, x)
> end
> 
> julia> @benchmark rand!(A)
> BenchmarkTools.Trial: 665 samples with 1 evaluation.
> Range (min … max): 7.472 ms … 9.067 ms ┊ GC (min … max): 0.00% … 0.00%
> Time (median): 7.492 ms ┊ GC (median): 0.00%
> Time (mean ± σ): 7.523 ms ± 140.133 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
> 
> ▇█▇▄▁
> █████▇█▇▆▆▄▅▆▅▅▁▅▅▅▄▁▄▆▁▁▄▁▅▁▄▄▄▄▁▄▁▁▁▁▄▄▄▄▁▅▁▁▁▁▁▁▁▁▁▄▁▄▄▅ ▇
> 7.47 ms Histogram: log(frequency) by time 8.14 ms <
> 
> Memory estimate: 0 bytes, allocs estimate: 0.
> 
> ```

### status quo:

```Julia-repl
julia> A = Array{Float32}(undef, 2^24);

julia> @benchmark rand!(A)
BenchmarkTools.Trial: 949 samples with 1 evaluation.
 Range (min … max): 5.200 ms … 6.274 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 5.240 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 5.268 ms ± 83.680 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

     ▃█▅▃▁▂
  ▄▅▅███████▆▇▆▆▆▇▅▆▆▆▄▆▅▇▆▆▅▄▄▅▁▅▄▄▅▄▁▄▄▁▁▄▄▄▁▄▁▁▁▁▄▁▄▁▁▄▁▄ ▇
  5.2 ms Histogram: log(frequency) by time 5.7 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

### After the above patch:

```Julia-repl
julia> @benchmark rand!(A)
BenchmarkTools.Trial: 665 samples with 1 evaluation.
 Range (min … max): 7.472 ms … 9.067 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 7.492 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 7.523 ms ± 140.133 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▇█▇▄▁
  █████▇█▇▆▆▄▅▆▅▅▁▅▅▅▄▁▄▆▁▁▄▁▅▁▄▄▄▄▁▄▁▁▁▁▄▄▄▄▁▅▁▁▁▁▁▁▁▁▁▄▁▄▄▅ ▇
  7.47 ms Histogram: log(frequency) by time 8.14 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

---

_[View the full topic](https://discourse.julialang.org/t/output-distribution-of-rand-float32-and-rand-float64-thread-2/105184)._