fork download
copy 
  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))
  101.  
Success #stdin #stdout 2.41s 37872KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
park_miller_bits  95/100
randu_bits        61/100
java_bits         89/100
python_bits       91/100
ansic_bits        95/100
https://ideone.com/2P7Xad
language:
Python 3 nbc (python 3.7.3)
created:
6 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!