diff --git a/Project.toml b/Project.toml
index 9841d72e..c021e5d0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -33,6 +33,7 @@ MvNormalCDF = "0.2, 0.3"
 Random = "1.6"
 Roots = "1, 2"
 SpecialFunctions = "2"
+StableRNGs = "1"
 StatsBase = "0.33, 0.34"
 TaylorSeries = "0.12, 0.13, 0.14, 0.15"
 Test = "1.6"
@@ -45,8 +46,9 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Test", "TestItemRunner", "InteractiveUtils", "LinearAlgebra", "HypothesisTests", "Aqua"]
+test = ["Test", "TestItemRunner", "InteractiveUtils", "LinearAlgebra", "HypothesisTests", "Aqua", "StableRNGs"]
diff --git a/test/archimedean_tests.jl b/test/archimedean_tests.jl
index cd5a1edd..f142d7f6 100644
--- a/test/archimedean_tests.jl
+++ b/test/archimedean_tests.jl
@@ -1,7 +1,8 @@
 
 @testitem "Test of τ ∘ τ_inv bijection" begin
-    using Random 
-    Random.seed!(123)
+    using Random
+    using StableRNGs
+    rng = StableRNG(123)
     taus = [0.0, 0.1, 0.5, 0.9, 1.0]
 
     for T in (
@@ -22,64 +23,76 @@ end
 
 
 @testitem "AMHCopula - Test sampling all cases" begin
+    using StableRNGs
+    rng = StableRNG(123)
     for d in 2:10
-        for θ ∈ [-1.0,-rand(),0.0,rand()]
-            rand(AMHCopula(d,θ),10)
+        for θ ∈ [-1.0,-rand(rng),0.0,rand(rng)]
+            rand(rng,AMHCopula(d,θ),10)
         end
     end
     @test true
 end
 @testitem "ClaytonCopula - Test sampling all cases" begin
-    using Random
-    Random.seed!(123) # to avoid weird bug on 1.6
-    rand(ClaytonCopula(2,-1),10)
+    using StableRNGs
+    rng = StableRNG(123)
+    rand(rng,ClaytonCopula(2,-1),10)
     for d in 2:10
-        for θ ∈ [-1/(d-1) * rand(),0.0,-log(rand()), Inf]
-            rand(ClaytonCopula(d,θ),10)
+        for θ ∈ [-1/(d-1) * rand(rng),0.0,-log(rand(rng)), Inf]
+            rand(rng,ClaytonCopula(d,θ),10)
         end
     end
     @test true
 end
 @testitem "FrankCopula - Test sampling all cases" begin
-    rand(FrankCopula(2,-Inf),10)
-    rand(FrankCopula(2,log(rand())),10)
+    using StableRNGs
+    rng = StableRNG(123)
+    rand(rng,FrankCopula(2,-Inf),10)
+    rand(rng,FrankCopula(2,log(rand(rng))),10)
     for d in 2:10
-        for θ ∈ [0.0,rand(),1.0,-log(rand()), Inf]
-            rand(FrankCopula(d,θ),10)
+        for θ ∈ [0.0,rand(rng),1.0,-log(rand(rng)), Inf]
+            rand(rng,FrankCopula(d,θ),10)
         end
     end
     @test true
 end
 @testitem "GumbelCopula - Test sampling all cases" begin
+    using StableRNGs
+    rng = StableRNG(123)
     for d in 2:10
-        for θ ∈ [1.0,1-log(rand()), Inf]
-            rand(GumbelCopula(d,θ),10)
+        for θ ∈ [1.0,1-log(rand(rng)), Inf]
+            rand(rng,GumbelCopula(d,θ),10)
         end
     end
     @test true
 end
 @testitem "JoeCopula - Test sampling all cases" begin
+    using StableRNGs
+    rng = StableRNG(123)
     for d in 2:10
-        for θ ∈ [1.0,1-log(rand()), Inf]
-            rand(JoeCopula(d,θ),10)
+        for θ ∈ [1.0,1-log(rand(rng)), Inf]
+            rand(rng,JoeCopula(d,θ),10)
         end
     end
     @test true
 end
 
 @testitem "GumbelBarnettCopula - Test sampling all cases" begin
+    using StableRNGs
+    rng = StableRNG(123)
     for d in 2:10
-        for θ ∈ [0.0,rand(),1.0]
-            rand(GumbelBarnettCopula(d,θ),10)
+        for θ ∈ [0.0,rand(rng),1.0]
+            rand(rng,GumbelBarnettCopula(d,θ),10)
         end
     end
     @test true
 end
 
 @testitem "InvGaussianCopula - Test sampling all cases" begin
+    using StableRNGs
+    rng = StableRNG(123)
     for d in 2:10
-        for θ ∈ [rand(),1.0, Inf]
-            rand(InvGaussianCopula(d,θ),10)
+        for θ ∈ [rand(rng),1.0, Inf]
+            rand(rng,InvGaussianCopula(d,θ),10)
         end
     end
     @test true
diff --git a/test/biv_Plackett.jl b/test/biv_Plackett.jl
index a6f6657e..7acf281c 100644
--- a/test/biv_Plackett.jl
+++ b/test/biv_Plackett.jl
@@ -60,9 +60,10 @@ end
 end
 
 @testitem "PlackettCopula Sampling" begin
-    using Random
+    using StableRNGs
+    rng = StableRNG(123)
     n_samples = 100
     C = PlackettCopula(0.8)
-    samples = rand(C, n_samples)
+    samples = rand(rng,C, n_samples)
     @test size(samples) == (2, n_samples)
 end
diff --git a/test/biv_clayton.jl b/test/biv_clayton.jl
index 9c9c31bd..7ae24736 100644
--- a/test/biv_clayton.jl
+++ b/test/biv_clayton.jl
@@ -27,15 +27,17 @@ end
 
 @testitem "sample" begin
     using StatsBase, Random
+    using StableRNGs
+    rng = StableRNG(123)
     n = 10^5
     θ = [-0.5, 2, 10]
     
     for ϑ in θ
         c = ClaytonCopula(2,ϑ)
         # if ϑ < 0
-        #     @test_broken rand(c,n)
+        #     @test_broken rand(rng,c,n)
         # else
-            u = rand(c, n)
+            u = rand(rng,c, n)
             @test corkendall(u') ≈ [1.0 Copulas.τ(c); Copulas.τ(c) 1.0] atol = 0.01
         # end
     end
diff --git a/test/coherence_archimedeans.jl b/test/coherence_archimedeans.jl
index 0e9c49a6..449d0f5e 100644
--- a/test/coherence_archimedeans.jl
+++ b/test/coherence_archimedeans.jl
@@ -1,7 +1,8 @@
 
 @testitem "Test of coherence for archimedeans" begin
     using HypothesisTests, Distributions, Random
-    Random.seed!(123)
+    using StableRNGs
+    rng = StableRNG(123)
     cops = (
         # true represent the fact that cdf(williamson_dist(C),x) is defined or not. 
         (AMHCopula(3,0.6), true),
@@ -24,17 +25,17 @@
         (InvGaussianCopula(3,8),true)
     )
     n = 1000
-    spl = rand(n)
-    spl2 = rand(n)
+    spl = rand(rng,n)
+    spl2 = rand(rng,n)
     for (C,will_dist_has_a_cdf_implemented) in cops
-        spl .= dropdims(sum(Copulas.ϕ⁻¹.(Ref(C),rand(C,n)),dims=1),dims=1)
+        spl .= dropdims(sum(Copulas.ϕ⁻¹.(Ref(C),rand(rng,C,n)),dims=1),dims=1)
         will_dist = Copulas.williamson_dist(C)
         if will_dist_has_a_cdf_implemented
             pval = pvalue(ExactOneSampleKSTest(spl, will_dist),tail=:right)
             @test pval > 0.01
         end
         # even without a cdf we can still test approximately:
-        spl2 .= rand(will_dist,n)
+        spl2 .= rand(rng,will_dist,n)
         pval2 = pvalue(ApproximateTwoSampleKSTest(spl,spl2),tail=:right)
         @test pval2 > 0.01
     end
diff --git a/test/fgmCopula_test.jl b/test/fgmCopula_test.jl
index a2a9fb87..dc6fbf91 100644
--- a/test/fgmCopula_test.jl
+++ b/test/fgmCopula_test.jl
@@ -20,11 +20,13 @@ end
 end
 
 @testset "FGMCopula PDF" begin
+    using StableRNGs
+    rng = StableRNG(123)
     examples = [
         ([0.1, 0.2, 0.3], [1.308876232, 1e-4], [0.1,0.2,0.5,0.4]),
         ([0.5, 0.4, 0.3], [1.024123232, 1e-4], [0.3,0.3,0.3,0.3]),
         ([0.1, 0.1], [0.01, 1e-4], 0.0),
-        ([0.5, 0.4], [1, 1e-4], rand()),
+        ([0.5, 0.4], [1, 1e-4], rand(rng)),
     ]
     
     for (u, expected) in examples
@@ -34,8 +36,10 @@ end
 end
 
 @testitem "FGMCopula Sampling" begin
+    using StableRNGs
+    rng = StableRNG(123)
     n_samples = 100
     F = FGMCopula(3,[0.1,0.2,0.3,0.4])
-    samples = rand(F, n_samples)
+    samples = rand(rng,F, n_samples)
     @test size(samples) == (3, n_samples)
 end
\ No newline at end of file
diff --git a/test/margins_uniformity.jl b/test/margins_uniformity.jl
index 5eb028d5..ea5bc7bc 100644
--- a/test/margins_uniformity.jl
+++ b/test/margins_uniformity.jl
@@ -1,6 +1,7 @@
 @testitem "Test samples have uniform maginals in [0,1]" begin
     using HypothesisTests, Distributions, Random
-    Random.seed!(1234)
+    using StableRNGs
+    rng = StableRNG(123)
     cops = (
         # true represent the fact that cdf(williamson_dist(C),x) is defined or not. 
         AMHCopula(3,0.6),
@@ -32,7 +33,7 @@
         nfail = 0
         d = length(C)
         @show C
-        spl = rand(C,n)
+        spl = rand(rng,C,n)
         @assert all(0 <= x <= 1 for x in spl)
         for i in 1:d
             @test pvalue(ApproximateOneSampleKSTest(spl[i,:], U),tail=:right) > 0.01 # quite weak but enough at these samples sizes to detect really bad behaviors.
diff --git a/test/old_precompile_workflow.jl b/test/old_precompile_workflow.jl
index 48e18d04..c58805fb 100644
--- a/test/old_precompile_workflow.jl
+++ b/test/old_precompile_workflow.jl
@@ -1,6 +1,8 @@
 
 @testitem "standard functionality test" begin
     using Random, Distributions
+    using StableRNGs
+    rng = StableRNG(123)
     biv_cops = [
         GaussianCopula([1 0.7; 0.7 1]),
         TCopula(2,[1 0.7; 0.7 1]),
@@ -12,11 +14,11 @@
         AMHCopula(2,0.7)
     ]
     for C in biv_cops
-        u = Random.rand(C,10)
+        u = Random.rand(rng,C,10)
         pdf(C,[0.5,0.5])
         cdf(C,[0.5,0.5])
         D = SklarDist(C,[Gamma(1,1),Normal(1,1)])
-        u = Random.rand(D,10)
+        u = Random.rand(rng,D,10)
         pdf(D,[0.5,0])
         cdf(D,[0.5,0])
     end
diff --git a/test/some_tests.jl b/test/some_tests.jl
index ae0f731b..d04b8308 100644
--- a/test/some_tests.jl
+++ b/test/some_tests.jl
@@ -2,8 +2,10 @@
     
     using Distributions
     using Random
+    using StableRNGs
+    rng = StableRNG(123)
     MyD = SklarDist(ClaytonCopula(3,7),(LogNormal(),Pareto(),Beta()))
-    u = rand(MyD,10000)
+    u = rand(rng,MyD,10000)
     rand!(MyD,u)
     fit(SklarDist{ClaytonCopula,Tuple{LogNormal,Pareto,Beta}},u)
     fit(SklarDist{GaussianCopula,Tuple{LogNormal,Pareto,Beta}},u)
@@ -24,7 +26,7 @@ end
 #     M1 = Beta(2,3)
 #     M2 = LogNormal(2,3)
 #     D = SklarDist(C,(M1,M2))
-#     X = rand(D,1000)
+#     X = rand(rng,D,1000)
 #     loglikelihood(D,X)
 #     fit(SklarDist{TCopula,Tuple{Beta,LogNormal}},X) # should give a very high \nu for the student copula. 
 # end
@@ -117,6 +119,8 @@ end
 @testitem "working measure" begin 
     using Distributions
     using Random
+    using StableRNGs
+    rng = StableRNG(123)
     
     for C in (ClaytonCopula(4,7.0),GumbelCopula(2, 1.2))
         d = length(C)
@@ -126,8 +130,8 @@ end
         @test Copulas.measure(C,u,v) >= 0
         
         for i in 1:d
-            u[i] = rand()
-            v[i] = u[i] + rand()*(1-u[i])
+            u[i] = rand(rng)
+            v[i] = u[i] + rand(rng)*(1-u[i])
         end
         @test Copulas.measure(C,u,v) >= 0
     end
@@ -140,8 +144,8 @@ end
         @test_broken Copulas.measure(C,u,v) >= 0
         
         for i in 1:d
-            u[i] = rand()
-            v[i] = u[i] + rand()*(1-u[i])
+            u[i] = rand(rng)
+            v[i] = u[i] + rand(rng)*(1-u[i])
         end
         @test_broken Copulas.measure(C,u,v) >= 0
     end
diff --git a/test/survivaltesting.jl b/test/survivaltesting.jl
index 461be8f3..961787e3 100644
--- a/test/survivaltesting.jl
+++ b/test/survivaltesting.jl
@@ -1,11 +1,13 @@
 @testitem "Testing survival stuff" begin
     using Distributions
+    using StableRNGs
+    rng = StableRNG(123)
     C = ClaytonCopula(2,3.0) # bivariate clayton with theta = 3.0
     C90 = SurvivalCopula(C,(1,)) # flips the first dimension
     C270 = SurvivalCopula(C,(2,)) # flips only the second dimension. 
     C180 = SurvivalCopula(C,(1,2)) # flips both dimensions.
 
-    u1,u2 = rand(2)
+    u1,u2 = rand(rng,2)
     c = cdf(C,[u1,u2])
     p = pdf(C,[u1,u2])
     @test cdf(C90,[1-u1,u2]) == c
diff --git a/test/williamson_test.jl b/test/williamson_test.jl
index f1ef321e..fdf1324b 100644
--- a/test/williamson_test.jl
+++ b/test/williamson_test.jl
@@ -1,7 +1,8 @@
 
 @testitem "williamson test" begin
     using Distributions, Random
-    Random.seed!(123)
+    using StableRNGs
+    rng = StableRNG(123)
     taus = [0.0, 0.1, 0.5, 0.9, 1.0]
 
     ϕ_clayton(x, θ) = max((1 + θ * x),zero(x))^(-1/θ)
@@ -13,6 +14,6 @@
         WilliamsonCopula(x -> ϕ_clayton(x,-0.3),2)
     )
     for C in Cops
-        x = rand(C,1000)
+        x = rand(rng,C,1000)
     end
 end