Haskell function for looking at next prime equivalence class bias, as described in https://arxiv.org/abs/1603.03720

next_prime_bias.hs

{--
We provide a function 'nextPrimeBias' that when called as
  nextPrimeBias a b n r
returns the probability that, among the first r primes, a prime that is equivalent to a mod n is followed by a prime that is equivalent to b mod n.

We see that these probabilities are not evenly distributed among equivalence classes, as described in https://arxiv.org/abs/1603.03720 :
nextPrimeBias 7 1 10 100000 = 0.25736558
nextPrimeBias 7 3 10 100000 = 0.27695382
nextPrimeBias 7 7 10 100000 = 0.144993
nextPrimeBias 7 9 10 100000 = 0.3206876
--}


-- get primes from optinised prime sieve
import Data.Numbers.Primes

-- reduce primes mod n
primesMod :: Integer -> [Integer]
primesMod n = map (\p -> mod p n) primes

-- quick 'n' dirty vectorisation
(+#) :: (Integer, Integer) -> (Integer, Integer) -> (Integer, Integer)
(a,b) +# (c,d) = (a+c,b+d)

-- main function
nextPrimeBias :: Integer -> Integer -> Integer -> Int -> Float
nextPrimeBias a b n r = ratio $ counts a b (take r $ primesModPairs)
  where
    primesModPairs = zip (primesMod n) $ tail (primesMod n)
    counts a b [] = (0,0)
    counts a b ((x,y):xys) | x /= a    = counts a b xys
                           | y == b    = (1,1) +# counts a b xys
                           | otherwise = (0,1) +# counts a b xys
    ratio pair = (fromIntegral $ fst pair) / (fromIntegral $ snd pair)