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:
- Repeatability (given an initial state, or seed, the output will be the same)
- Number of values (a period is the maximum number of states a generator can make, a bigger period means more random values before repeating)
- 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 = (aXn + 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 = (aXn + cn-1) mod b, cn = (aXn + 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)
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.
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 (
mt, the ‘cache’
mti, ‘cache’ pointer (which cache value to get)
N, the size of the ‘cache’ (
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 y = (self.mt[self.N-1] & self.UPPER_MASK) | (self.mt & self.LOWER_MASK) 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()