Adds test that shows thread-safety improvement for rhyper

test.jl

@@ -12,3 +12,48 @@ unsafe_store!(cglobal((:exp_rand_ptr, libRmath), Ptr{Cvoid}),
     @test ccall((:dbeta, libRmath), Float64, (Float64, Float64, Float64, Int32), 0.5, 0.1, 5.0, 0) ≈ 0.014267678091051986
     @test 0 <= ccall((:rbeta, libRmath), Float64, (Float64, Float64), 0.1, 5.0) <= 1.0
 end
+
+@testset "rhyper" begin
+    # double rhyper(double nn1in, double nn2in, double kkin)
+    Nred = 30.0
+    Nblue = 40.0
+    Npulled = 5.0
+
+    hyper_samples = [
+        ccall((:rhyper, libRmath), Float64, (Float64, Float64, Float64), Nred, Nblue, Npulled)
+        for _ in 1:1_000_000
+    ]
+    expected_mean = Npulled * Nred / (Nred + Nblue)
+    sample_mean = sum(hyper_samples) / length(hyper_samples)
+    @test sample_mean ≈ expected_mean rtol = 0.01
+
+    N = (Nred + Nblue)
+    expected_variance = Npulled * Nred * (N - Nred) * (N - Npulled) / (N * N * (N - 1))
+    sample_variance = 1 / (length(hyper_samples)) * sum((hyper_samples .- sample_mean) .^ 2)
+    @test sample_variance ≈ expected_variance rtol = 0.01
+end
+
+using Distributions
+function sample_KkC(n; N, Q)
+    total_errors = Distributions.Binomial(N, Q)
+    K = rand(total_errors)
+    k = ccall(
+        (:rhyper, libRmath), Float64, (Float64, Float64, Float64),
+        K, N-K, n
+    )
+    return k
+end
+
+@testset "fulll" begin
+    function f(Q)
+        objective(n) = [sample_KkC(n; N = 819_200, Q) for _ = 1:100]
+        vals = [10, 100]
+        objective.(vals)
+    end
+
+    Qs = [0.05, 0.055, 0.1, 0.2, 0.3]
+
+    Threads.@threads for i in eachindex(Qs)
+        f(Qs[i])
+    end
+end