Random Number Generators

Random number generators are used throughout programming and game design. So I wanted to know how does a random number generator work. I looked into quite a few generator algorithms and wrote my own generators, in python (C# too but that’s on github only), for what I considered to be the main ones. One thing I quickly noticed from my research is these aren’t true random number generators (TRNG) but pseudo-random number generators (PRNG), the numbers generated only approximate randomness. There are also random number generators used for cryptography (CPRNG) that I didn’t look into. For more information about the difference of TRNG and PRNG you can read this article.

The random number algorithms I examined were:

• Linear Congruential Generator (LCG) – one of the oldest and best-known PRNG
• Blum Blum Shub Generator (BBS) – probably the coolest name for an algorithm
• Multiply with Carry (MWC)
• Xorshift (XOR)
• Mersenne Twister (MT) – most widely used general-purpose PRNG

Keep in mind these are only a few of the random number generator algorithms out there, and some of the listed algorithms have variations or improvements. Some of the varied algorithms from above:

• LCG -> Combined Linear Congruential Generator
• MWC-> Complementary Multiply with Carry (informally known as the Mother of All PRNGs)
• Xorshift-> Xorshift*, Xorshift+, Xoroshiro128
• MT-> WELL215

Some other random number generator algorithms:

• Middle Square Method
• Fibonacci-Lagged generator
• PCG Family

The goal of a good PRNG is:

1. Repeatability (given an initial state, or seed, the output will be the same)
2. Number of values (a period is the maximum number of states a generator can make, a bigger period means more random values before repeating)
3. Speed (fast to compute)

Choosing a good seed, initial state, is very important in random number algorithms. Some seeds will give very short periods. Also choosing the same seed will give the same sequence of numbers every time. Most often the seed given is the number of seconds since the epoch, 1 January 1970. There are some interesting alternatives to generating a random seed besides using time or user entered seed; using the PID (process id, good idea for multithreading) or user input (PuTTY has the user draw randomly/move the mouse around to generate a SSH key or entropy -randomness generated by the computer based on user input).

Now the code. I will show the functions and if you want to see them all together I’ve included them on my github. If you want to return a number in a range [a,b], inclusive, you can write a function like:

```'''
mini-minimum value
maxi-maximum value

next_PRN (next pseudo random number)-use whichever algorithm you want
'''
def next_RN_range(mini, maxi):
return mini + (int)(next_PRN() * ((maxi - mini) + 1))
```

To return a decimal in all of my code I divide the random number by another, usually the modulus or highest number.

Linear Congruential Generator

One of the simplest generators to understand because it has a very simple equation.

Xn+1 = (aX+ c) mod m

a, 0 < a < m, the multiplier (arbitrary, Lehmer used 23)

c,  0 <= c < m, the increment (arbitrary 1, 11, 12345)

m, 0 < m, the modulus (usually based on a power of 2, 232, 264, 231-1)

X0, 0<= X0 <= m, the seed (avoid 0, time since the epoch is commonly used)

```class RandomGenerators:
def __init__(self, seed = None):
'''
If you want to understand why having a bigger period is important, use the following
numbers which produce a small period:
self.lcg_mod = 10
self.lcg_a = 2
self.lcg_c = 3
'''
self.lcg_mod = 2**31
self.lcg_a = 1103515245
self.lcg_c = 12345
# seed
if seed == None:
self.lcg_prev = int(time.time()) # seed set to number of seconds since epoch
else:
self.lcg_prev = seed

def next_LCG(self):
rand_num = (self.lcg_prev * self.lcg_a + self.lcg_c) % self.lcg_mod
self.lcg_prev = rand_num # update the state
return self.lcg_prev / float(self.lcg_mod)
```

Multiply with Carry

This is similar to LCG but the c is also calculated based on the previous c, so the last c (in addition to the seed/last-random-number) has to be kept track of from one state to the next.

Xn+1 = (aX+ cn-1) mod b, cn = (aX+ cn-1) / b

a, multiplier (arbitrary)

b, base (usually half of max double/int)

X0, 0<= X0 <= m, the seed (avoid 0, time since the epoch is commonly used)

c0, the carry (arbitrary)

```class RandomGenerators:
def __init__(self, seed = None, carry = None):
self.mwc_b = 2**32      # base, usually half of double/int size
self.mwc_a = 1103515245 # multiplier
# seed
if seed == None:
self.mwc_prev = int(time.time()) # seed set to number of seconds since epoch
else:
self.mwc_prev = seed
#carry
if carry == None:
self.mwc_c = 3 # arbitrary carry
else:
self.mwc_c = carry

def next_MWC(self):
t = self.mwc_a * self.mwc_prev + self.mwc_c
rand_num = t % self.mwc_b
self.mwc_prev = rand_num
self.mwc_c = int(t / self.mwc_b)
return rand_num / float(self.mwc_b)```

Blum Blum Shub

This is another fairly straightforward algorithm that has a max period of length M.

Xn+1 = Xn2 mod M

M, modulus (product of two large primes p and q; both of which should be congruent to 3, ie… (p or q) % 4 == 3)

X0, the seed (should be coprime with M…have it be a prime for simplicity, can be another number)

```class RandomGenerators:
def __init__(self):
p = 492876847
q = 715225739
self.bbs_mod = p * q
self.bbs_prev = 920419823

def next_BlumBlumShub(self):
rand_num = (self.bbs_prev**2) % self.bbs_mod
self.bbs_prev = rand_num
return rand_num / float(self.bbs_mod)```

Xorshift

This PRNG doesn’t use an equation, instead, it uses bit operations. It xor, exclusive-or, the current random number with the current random number bit shifted left some number. Then it does another xor but with a right bit shift. And then one final time, again left bit shifted. In my example I use 21, 35, and 4 but other numbers work (some other sets that have good results are (12, 25, 27), (13, 7, 17), (13, 17, 5)).

X0, the seed

```class RandomGenerators:
def __init__(self, seed = None):
# just need a seed
if seed == None:
self.xor_prev = int(time.time())
else:
self.xor_prev = seed

def next_xorshift(self):
rand_num = self.xor_prev
rand_num ^= (rand_num << 21)
rand_num ^= (rand_num >> 35)
rand_num ^= (rand_num << 4)
self.xor_prev = rand_num
return (rand_num & 0xFFFFFFFF) / float(0xFFFFFFFF)```

In the return: `& 0xFFFFFFFF` truncates, or grabs the 8 rightmost bits. `^` stands for xor, or exclusive or (1 if only one of the two compared bits is 1, not if both are, and not if neither are). `<<` is left shift and `>>` is right shift.

Mersenne Twister

This was the most complicated PRNG that I looked at. If you break it down it reminds me a bit of the XOR shift. I’ll add comments in the code to try to explain what’s going on but mostly its bit operations. This can be seeded with an array ( `init_by_array()`) or a single number (`init_genrand()`).

`mt`, the ‘cache’

`mti`, ‘cache’ pointer (which cache value to get)

`N`, the size of the ‘cache’ (`mt`)

A common snippet is `&= 0xFFFFFFFF` which will get the 8 rightmost bits. My random numbers do not match those produced by the C code because python’s hexadecimal doesn’t do negatives the same way. I don’t know how this affects the period or randomness of the generator, though.

```class RandomGenerators:
# period parameters
N = 624
M = 394
MATRIX_A = 0x9908b0df   # constant vector a
UPPER_MASK = 0x80000000 # most significant w-r bits
# this will get the most significant bit (MSB), leftmost
LOWER_MASK = 0x7fffffff # least significant r bits
# this will get the least significant bits, rightmost

def __init__(self):
self.mt = [None]
self.mti = self.N + 1

'''
initialize by an array with array-length
init_key is the array for initializing keys
key_length is its length
'''
def init_by_array(self, init_key, key_length):
# fill array with arbitrary seed
self.init_genrand(19650218)
i,j = 1,0
# k gets the higher value: N or key_length
k = self.N if self.N > key_length else key_length
while k > 0:
self.mt[i] = (self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 30)) * 1664525)) + init_key[j] + j # non linear
self.mt[i] &= 0xFFFFFFFF
i+=1
j+=1
if i >= self.N:
self.mt = self.mt[self.N-1]
i = 1
if j >= key_length:
j = 0
k-=1

k = self.N-1
while k > 0:
self.mt[i] = (self.mt[i] ^ ((self.mt[i-1] ^ (self.mt[i-1] >> 30)) * 1566083941)) - i # non linear
self.mt[i] &= 0xffffffff

i += 1
if i >= self.N:
self.mt = self.mt[self.N-1]
i=1
k-=1

self.mt = 0x80000000 # MSB is 1; assuring non-zero initial array

# Adding this makes mt have the same values the C does at this point
#for i,num in enumerate(self.mt):
#if num>=2147483648: # (2**32)/2
#self.mt[i]-=4294967296 # 2**32''''

# initializes mt[N] with a seed
def init_genrand(self,s):
self.mt = s & 0xffffffff
self.mti = 1
while self.mti < self.N:
self.mt.append(1812433253 * (self.mt[self.mti-1] ^ (self.mt[self.mti - 1] >> 30)) + self.mti)
self.mt[self.mti] &= 0xffffffff
self.mti += 1

# generates a random number on [0,0xffffffff]-interval
def genrand_int32(self):
y = 0
mag01=[0x0, self.MATRIX_A]
if self.mti >= self.N: # generate N words at one time
kk = 0
if self.mti == self.N + 1: # if init_genrand() has not been called,
self.init_genrand(5489) # a default initial seed is used

while kk < self.N - self.M:
y = (self.mt[kk] & self.UPPER_MASK) | (self.mt[kk+1] & self.LOWER_MASK)                 self.mt[kk] = self.mt[kk+self.M] ^ (y >> 1) ^ mag01[y & 0x1]
kk += 1

while kk < self.N - 1:
y = (self.mt[kk] & self.UPPER_MASK) | (self.mt[kk+1] & self.LOWER_MASK)                 self.mt[kk] = self.mt[kk+(self.M-self.N)] ^ (y >> 1) ^ mag01[y & 0x1]
kk += 1

self.mt[self.N-1] = self.mt[self.M-1] ^ (y >> 1) ^ mag01[y & 0x1]

self.mti = 0

y = self.mt[self.mti]
self.mti += 1

# Tempering
y ^= (y >> 11)
y ^= (y << 7) & 0x9d2c5680
y ^= (y << 15) & 0xefc60000
y ^= (y >> 18)

return y

# generates a random number on [0,0x7fffffff]-interval
def genrand_int31(self):
return (self.genrand_int32() >> 1)

# generates a random number on [0,1]-real-interval
def genrand_real1(self):
return self.genrand_int32() * (1.0 / 4294967295.0)
# divided by 2^32-1

# generates a random number on [0,1)-real-interval
def genrand_real2(self):
return self.genrand_int32() * (1.0 / 4294967296.0)
# divided by 2^32

# generates a random number on (0,1)-real-interval
def genrand_real3(self):
return ((self.genrand_int32()) + 0.5) * (1.0 / 4294967296.0)
# divided by 2^32

# generates a random number on [0,1) with 53-bit resolution
def genrand_res53(self):
a, b = self.genrand_int32()>>5, self.genrand_int32()>>6
return(a * 67108864.0 + b) * (1.0 / 9007199254740992.0)

r = RandomGenerators()
init = [0x123, 0x234, 0x345, 0x456]
length = 4
r.init_by_array(init, length)
for i in range(10):
print(str(r.genrand_int31())+' ', end="")
if i%5==4:
print()
print()
for i in range(10):
print(str(r.genrand_real2())+' ', end="")
if i%5==4:
print()```