Random numbers generator

A long time ago, in the early eighties, I had a one of the legendary home computers. That computer was my introduction to game programming. Despite the games I created, I was always fascinated by one specific game. The time that has passed since the last time I played it was enough to make me forget its name. But I never forgot what excited me the most.

It was a shoot them up game where you launched a space ship from a space station and your goal was to shoot the enemy ships before they destroyed the base. The most interesting part for me though was that the stars were always at the same position. Their number was very big to hold all the data in memory considering that it was only 32 kilobytes for the language interpreter, program, data AND video ram. They seemed to be spread realistically random in the sky. So what was the trick?

Years later as I was doing some research work for a cryptography project I was working on I found something that reminded me of the stars distribution problem. It was a random number generator. It was a small and efficient algorithm that was generating sequences of numbers in an unpredictable manner and good distribution over a given domain.

Computers cannot produce real random numbers. The very nature of the device, it is actually a deterministic device, prohibits this. They can only generate long sequences of seemingly random numbers. The statistical distribution of these numbers should be rather satisfactory and any number within the given domain should be equally possible to occur. Another factor we must take into account is unpredictability. This means that although the sequence (1, 2, 3 ... 100) is statistically perfect it is very predictable and cannot be considered random.

Linear Congruential Method

The most well-known method for generating random numbers, which has been used almost exclusively since it was introduced by D. Lehmer in 1951, is the so-called linear congruential method. It is a very straight forward method that does not take up many resources. Let us assume that seed contains some arbitrary initial value. We will use the following code to generate 'random' numbers.

static int a;
static int b;
static int m;
void init_random(int seed)
{
	a = seed;
}
int random()
{
	a = (a*b+1) % m;
	return a;
}

In the code fragment above we introduce the variables 'a', 'b' and m. The value of 'a' changes with every call to the random function, based on the constants 'b' and 'm'. The result is always between 0 and m-1.

Although this seems simple there are a couple of points we have to pay some attention to. First is the selection of 'seed', 'b' and 'm'. First 'm' should be large, close to the computer word size. It is also convenient to make it a power of 2. Second 'b' should be neither too big nor too small. A good choice is to make it one digit smaller than 'm'. Third you should make sure that it ends with the pattern 'x21' where 'x' is even. This has been analyzed by many mathematicians over the years.

Finally we have to consider the fact that computer calculations have a simple but quite significant characteristic that cannot be ignored: overflow. This will eventually occur in the multiplication of 'a' and 'b' and will lead into unexpected results. So we write the 'safe_multiply' function that breaks the multiplication into smaller parts less likely to overflow. The principle is elementary

pq = (10^4 * p1+p0) (10^4 * q1+q0) = 10^8 * p1*q1 + 10^4(p1*q0+p0*q1) + p0*q0

#include <stdio.h>
static int a = 0;
const int b = 32456821; // notice the ..x21 pattern
const int m = 100000000; // the maximum +1

void init_random(short seed)
{
	a = seed;
}

// multiply eliminating overflow
int safe_multiply(int p, int q)
{
	const int m1 = 10000;
	int p1,p0,q1,q0;
	int ret;
	p1 = p/m1; p0 = p%m1;
	q1 = q/m1; q0 = q%m1;

	ret = ((p0*q1+p1*q0)%m1)*m1+p0*q0;

	return ret;
}
int random()
{
	a = (safe_multiply(a,b)+1) % m;
	return a;
}

int main(int argc, char** argv)
{
	int i;
	init_random(12345);
	// generate 20 random numbers in the range 0-99
	for (i=0; i<20; i++)
		printf("%d -> %d\n", i, random()%100);
	return 0;
}