diff --git a/Project.toml b/Project.toml index ae56068..9a55e22 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,18 @@ -name = "StableRNGs" -uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" +name = "StableRNGs" +uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" authors = ["Rafael Fourquet "] -version = "1.0.1" +version = "1.0.2" [deps] Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] +Random = "<0.0.1, 1" +Test = "<0.0.1, 1" julia = "1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/src/StableRNGs.jl b/src/StableRNGs.jl index c661bc2..af62491 100644 --- a/src/StableRNGs.jl +++ b/src/StableRNGs.jl @@ -155,4 +155,29 @@ function Random.shuffle!(r::StableRNG, a::AbstractArray) return a end +# https://github.com/JuliaRandom/StableRNGs.jl/issues/20 +@noinline function Random.randn_unlikely(rng::StableRNG, idx, rabs, x) + @inbounds if idx == 0 + while true + xx = -Random.ziggurat_nor_inv_r*log(rand(rng)) + yy = -log(rand(rng)) + yy+yy > xx*xx && + return (rabs >> 8) % Bool ? -Random.ziggurat_nor_r-xx : Random.ziggurat_nor_r+xx + end + elseif (Random.fi[idx] - Random.fi[idx+1])*rand(rng) + Random.fi[idx+1] < exp(-0.5*x*x) + return x # return from the triangular area + else + return randn(rng) + end +end +@noinline function Random.randexp_unlikely(rng::StableRNG, idx, x) + @inbounds if idx == 0 + return Random.ziggurat_exp_r - log(rand(rng)) + elseif (Random.fe[idx] - Random.fe[idx+1])*rand(rng) + Random.fe[idx+1] < exp(-x) + return x # return from the triangular area + else + return Random.randexp(rng) + end +end + end # module diff --git a/test/runtests.jl b/test/runtests.jl index e437e65..f7a9165 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -140,3 +140,30 @@ end shuffle!(StableRNG(10), b) @test b == a_shuffled end + +# https://github.com/JuliaRandom/StableRNGs.jl/issues/20 +@testset "`randn` stability" begin + ref = [ + 0.5745734638645761, + 0.9050768627399978, + 0.7998353512850861, + 3.8845391427592286, + -0.9209167676765456, + -0.8486914352114853, + -1.187370886634173 + ] + @test randn(StableRNG(1_337), 10_000)[4214:4220] == ref +end + +@testset "`randexp` stability" begin + ref = [ + 1.7805339657229489, + 4.29332694381576, + 0.20777989530218552, + 8.196071589366719, + 7.551925528256079, + 1.3540162045313204, + 0.5239664874260928 + ] + @test randexp(StableRNG(1_337), 10_000)[1545:1551] == ref +end