diff --git a/src/fit.jl b/src/fit.jl
index 83e657c5..a132752e 100644
--- a/src/fit.jl
+++ b/src/fit.jl
@@ -560,21 +560,26 @@ end
 mutable struct RandomCoupels
     n::Int
     count::Int
+    iseed::Int
 end
 
-function Base.iterate(iter::RandomCoupels, state = 0)
-    if state == iter.count
+Base.length(iter::RandomCoupels) = iter.count
+
+function Base.iterate(iter::RandomCoupels, state = (0, MersenneTwister(iter.iseed)))
+    count,rng = state
+
+    if count == iter.count
         return nothing
     end
 
     # pick two random points
-    j = rand(1:iter.n)
+    j = rand(rng,1:iter.n)
     i = j
     while (i == j)
-        i = rand(1:iter.n)
+        i = rand(rng,1:iter.n)
     end
 
-    return ((i, j), state + 1)
+    return ((i, j), (count + 1, rng))
 end
 
 
@@ -652,6 +657,7 @@ this function used to be called lfit in fitlsn.f
 function fitlen(x::Tuple, d, weight, nsamp; kwargs...)
     # number of samples
     n = length(d)
+    iseed = n
 
     iter =
         if (nsamp == 0)
@@ -662,7 +668,7 @@ function fitlen(x::Tuple, d, weight, nsamp; kwargs...)
                 @warn "Strange to ask for more samples than available from data; will proceed"
             end
 
-            RandomCoupels(n, (nsamp * (nsamp - 1)) ÷ 2)
+            RandomCoupels(n, (nsamp * (nsamp - 1)) ÷ 2,iseed)
         end
 
     if (n > 10000) && (nsamp != 0)
@@ -1035,8 +1041,8 @@ function fithorzlen(
 
     weight = 1 ./ epsilon2
 
-    #Threads.@threads for k = 1:length(z)
-    for k = 1:length(z)
+    Threads.@threads for k = 1:length(z)
+    #for k = 1:length(z)
 
         sel = if length(x) == 3
             (abs.(x[3] .- z[k]) .< searchz)
diff --git a/test/test_fit.jl b/test/test_fit.jl
index f4edf04c..b1a1a753 100644
--- a/test/test_fit.jl
+++ b/test/test_fit.jl
@@ -1,5 +1,7 @@
 using Test
 using DelimitedFiles
+using Statistics
+using Random
 import DIVAnd
 
 # test data for basic statistics
@@ -167,3 +169,8 @@ fitlenz, dbinfo = @test_logs (:info, r".*at*") match_mode = :any DIVAnd.fitvertl
     epsilon2 = epsilon2,
 );
 @test median(fitlenz) ≈ lenz rtol = 0.5
+
+
+# RandomCoupels iterators
+
+@test collect( DIVAnd.RandomCoupels(1000,10,1234) ) == collect( DIVAnd.RandomCoupels(1000,10,1234) )