search

JuliaSIMD / VectorizedRNG.jl

Vectorized uniform and normal random samplers.
MIT License
33 stars 7 forks source link

Bug of randn for Float32? #19

Closed AStupidBear closed 2 years ago

AStupidBear commented 2 years ago 
julia> using Random, Statistics

julia> using VectorizedRNG

julia> rng = local_rng()

VectorizedRNG.Xoshift{2}(Ptr{UInt64} @0x00000000011c52c0)

julia> T, N = Float32, 100000
(Float32, 100000)

julia> std(randn(rng, T, N))
0.81086445f0

julia> T, N = Float64, 100000
(Float64, 100000)

julia> std(randn(rng, T, N))
0.9988777670204232

julia> versioninfo()
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) Platinum 8255C CPU @ 2.50GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, cascadelake)
Environment:
  JULIA_EDITOR = code
chriselrod commented 2 years ago

I can reproduce. My nlogn function is wrong for Float32. Replacing it with -sqrt(-log(floatbitmask(u2,T)-oneopenconst(T))) produces

julia> std(randn(rng, T, N))
1.0006113f0

julia> std(randn(rng, T, N))
1.0034106f0

julia> std(randn(rng, T, N))
0.99588895f0

julia> std(randn(rng, T, N))
1.000607f0

Might be a good idea to translate this algorithm to be architecture generic and use Float32 to avoid the division the SLEEFPirates uses: https://github.com/JuliaSIMD/VectorizationBase.jl/blob/7adea842ca9bbaed371222ac863f90ce302db6f6/src/special/log.jl#L215-L254 Can of course use a smaller polynomial for Float32.

AStupidBear commented 2 years ago

Could you please temorarily update the code? I have no idea how to patch it.

chriselrod commented 2 years ago

Could you please temorarily update the code? I have no idea how to patch it.

Sure, fixed by https://github.com/JuliaSIMD/VectorizedRNG.jl/commit/c3bafaf83f3e511081b9b4d51b98e950496cc2aa.