ploffie

prng tests

Jul 10th, 2017
119
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. from math import log, exp, sqrt, floor, ceil, log2
  2. import itertools as IT
  3. import random
  4. from collections import Counter
  5.  
  6. MOD = 2**31
  7.  
  8. def rint(m):
  9.     return random.randrange(0, m)
  10.  
  11. def ansic(seed=None):
  12.     M = 2**31
  13.     if seed is None:
  14.         seed = rint(0, M)
  15.     while True:
  16.         seed = (seed * 1103515245 + 12345) % M
  17.         yield seed >> 16
  18.  
  19. def randu(seed=None):
  20.     if seed is None:
  21.         seed = rint(MOD)
  22.     while True:
  23.         seed = (seed * 65539) % MOD
  24.         yield seed
  25.  
  26. def park_miller(seed=None):
  27.     if seed is None:
  28.         seed = rint(MOD)
  29.     while True:
  30.         seed = (seed * 16807) % (MOD-1)
  31.         yield seed
  32.  
  33. def java(seed=None):
  34.     F = 2**48
  35.     if seed is None:
  36.         seed = rint(F)
  37.     while True:
  38.         seed = (seed * 25214903917 + 11) % (F-1)
  39.         yield seed >> 16
  40.  
  41. def prng_to_bits(rng, nbits):
  42.     while True:
  43.         b = bin(next(rng))[2:]
  44.         yield from (int(bi) for bi in reversed(b))
  45.         yield from IT.repeat(0, nbits - len(b))
  46.  
  47. def create_bits(rng, nbits, m):
  48.     def bits():
  49.         yield from prng_to_bits(rng(rint(m)), nbits)
  50.     bits.__name__ = rng.__name__ + "_bits"
  51.     return bits
  52.  
  53. park_miller_bits = create_bits(park_miller, 31, MOD)
  54. randu_bits = create_bits(randu, 31, MOD)
  55. ansic_bits = create_bits(ansic, 15, MOD)
  56. java_bits = create_bits(java, 32, 2**48)
  57.  
  58. def python_bits():
  59.     while True:
  60.         yield random.getrandbits(1)
  61.  
  62. def count_bits(n):
  63.     return bin(n).count("1")
  64.  
  65. def ones_test(prng, n=200):
  66.     ones  = sum(count_bits(next(prng)) for _ in range(n))
  67.     expected = n / 2
  68.     stdev = sqrt(n / 4)
  69.     lo, hi = round(expected - 2 *stdev), round(expected + 2 * stdev)
  70.     if not lo < ones < hi:
  71.         # print("Ones not {} < {} < {} - {}".format(lo, ones, hi, prng.__name__))
  72.         return 0
  73.     return 1
  74.  
  75. def count_bit_runs(prng, n=1000):
  76.     return Counter(len(list(g)) for
  77.                    k, g in IT.groupby(IT.islice(prng, n)) if k == 1)
  78.  
  79. def runs_test(prng, nbits=100000):
  80.     runs = count_bit_runs(prng, nbits)
  81.     longest_run = max(runs.keys())
  82.     eruns, esdev =  -log(nbits/2)/log(0.5), 1 / log(2)
  83.     rlo, rhi = ceil(eruns - 2 * esdev), floor(eruns + 2 * esdev)
  84.     if not rlo <= longest_run <= rhi:
  85.         # print("Run not {} < {} < {} - {}".format(rlo, longest_run, rhi, prng.__name__))
  86.         return 0
  87.     return 1
  88.  
  89. def tests(prng, nbits=100000):
  90.     s1 = ones_test(prng, n=nbits)
  91.     s2 = runs_test(prng, nbits)
  92.     return s1 + s2
  93.  
  94. N = 10000
  95. ncases = 50
  96. for rng in (park_miller_bits, randu_bits, java_bits, python_bits, ansic_bits):
  97.     score = 0
  98.     for i in range(ncases):
  99.         score += tests(rng(), N)
  100.     print("{:16s}  {}/{}".format(rng.__name__, score, 2*ncases))
RAW Paste Data