Fix streams of randn and randexp on Julia >= 1.11 (#22)

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 <[email protected]>"]
-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"]

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

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