Page 1 of 1

Random Number Generation 102 Coding a Linear Congruential Generator Rate Topic: -----

#1 NickDMax  Icon User is offline

  • Can grep dead trees!
  • member icon

Reputation: 2250
  • View blog
  • Posts: 9,245
  • Joined: 18-February 07

Posted 21 February 2007 - 07:11 AM

Generating Pseudorandom Numbers in C/C++
Coding a Linear Congruential Generator

This tutorial shows how to generate random numbers using a Linear Congruential Sequences. Most PRNGs used in programming languages use an LCG so you might as well just use rand() to generate your random data. However, maybe you are just interested in how these things work. Maybe you don’t trust the programmers who wrote that library. Maybe you are just a plane old-fashioned do-it-yourselfer.

I. Theory

The heart of an LCG is the following formula

X(i+1) = (a * X(i) + c) mod M

where

M is the modulus. M > 0.
a is the multiplier, 0 <= a < M.
c is the increment, 0 <= c < M.
X(0) is the seed value, 0 <= X(0) < M.
i is the iterator. i < M

Suppose we were to use the following:
#include <iostream>
using namespace std;

int main()
{
	int M = 8;
	int a = 5;
	int c = 3;
	int X = 1;
	int i;
	for(i=0; i<8; i++)
	{
		X = (a * X + c) % M;
		cout << X << “ “;
	}
	return 0;
} 

Output: 0 3 2 5 4 7 6 1

This generates a sequence of numbers from 0 to 7 in a fairly scrambled way. If we were to extend the loop we would find that this sequence would repeat over and over. This sequence happens to have a period of M, but other choices for a, and c could have smaller periods (a=3, c=3 has a period of 4).

To make our LCG useful we will need a large period. To do this we have to choose appropriate values for M, a, and c.

To have a maximum period of M the following conditions must be met:

1. M and c are Relatively Prime so the gcd(M, c) = 1.
2. (a-1) is divisible by all prime factors of M.
3. (a-1) mod 4 = 0 if (M mod 4) = 0

The first condition is the hardest to check. This is a link to a good GCD algorithm written in C/C++, however, there is no need to add this to our code unless we want to generate c and M on the fly. For the most part once we select a good set of values for c and M we can leave them since there are some really bad choices for c and M.

The second condition can be difficult for large values of M since there is no fast way to factor a large number. However if we manufacture M from prime numbers then we can easily find a value for a. In point of fact for ease of calculation M is often chosen to be a power of 2 and so this condition becomes “a is an odd number greater then 1 and less then M.” Easy enough.

The last condition is the easiest to verify. We simply need to divide M by 4 and if it divides evenly then we need to see if (a-1) also divides evenly by 4. Which is exactly what the Mod (%) operator does for us.

II. Failings

Armed with knowledge lets make a sequence 256 elements long. So M=2^8. This makes choosing c simple as it can be any odd number less than M. Lets choose c=27. Now choosing a value for a is also easy, choose a multiple of 4 and add 1. a = 65.

#include <iostream>
using namespace std;

int main()
{
	int M = 256;
	int a = 65;
	int c = 27;
	int X = 1;
	int i;
	for(i=0; i<256; i++)
	{
		X=(a * X + c) % M;
		cout << X << " ";
	}
	return 0;
}


Some of the output:
92 119  82 237  72  99  62 217  52  79
 42 197  32  59  22 177  12  39   2 157
248  19 238 137 228 255 218 117 208 235

Now if you are astute you will notice that this would not make a very good random number generator. Do you see it? The sequence is even, odd, even, odd, etc.. This is not a very *random* sequence. If we were to look at the numbers in binary we would notice other patterns. This is due to our choices of a, c, and M. Most importantly, if we choose a power of 2 for M then the least significant digits will have an even period less then M.

Another statistical anomaly of LCGs is found when we use the LCG to generate triplets in a 3-dimensional space. We find will find the points arrange themselves into planes.

At this point you are probably not sold on LGCs. There are however a few tricks that can be used to ameliorate many of these shortcomings.

1. Chose our constant wisely! A list of values used in some classical random number generators. (Such as the one used in ASCII-C’s rand() function). Some are good, some are bad!

2. Only use a few of the most significant bits X. This hides all of the high frequency patterns found in the sequence. If we choose a large value for M then we can be assured that we have a period that is larger than our expected data set.

3. Use obfuscation methods such as choosing variable numbers of bits or re-seeding.

Even with all of these it should be noted that LCG’s are not fit for cryptographic use, nor should they be used in statistical tests (yes, that’s right, all of those “coin-flipping” programs you wrote were using an unfair coin).

III Implementation

Here is a basic LCG random number generator:

/**********************************************
/*   A Basic Linear Congruential Generator
/*  ---------------------------------------
/*  Program By: NickDMax for Dreamincode.net
/*  02/23/2007
/*
/*  You may use this code as you see fit
/**********************************************/
#include <iostream>
#include <ctime>
using namespace std;

//This class will give satisfy basic
// random number considerations for
// games an general apps.
// WARNING: Not for Cryprographic use
//   Not for statistical methods.
class BasicLCG {
	
	private:
		unsigned long iCurrent;
	public:
		BasicLCG();
		BasicLCG(unsigned long);
		void seed(unsigned long iSeed);
		unsigned long nextNumber(); //get the next random number
		unsigned short int nextInt();
		unsigned char nextChar();
		int nextBit();
		double nextDouble();
		int inRange(int min, int max);
};

//Just a little test code to print some numbers
int main()
{
	BasicLCG rng(time(NULL));
	int i;
	//Lets see some bits...
	for( i=1; i<81; i++) {
		cout << rng.nextBit() <<"\t";
	}
	cout <<"\n";
	for( i=1; i<41; i++) {
		cout << (int)rng.nextChar() <<"\t";
	}
	cout <<"\n";
	for( i=1; i<41; i++) {
		cout << rng.nextInt() <<"\t";
	}
	cout <<"\n";
	for( i=1; i<41; i++) {
		cout << rng.nextNumber() <<"\t";
	}
	cout <<"\n\n";
	for( i=1; i<41; i++) {
		cout << rng.nextDouble() <<"\t";
	}
	cout <<"\n\n";
	for( i=1; i<41; i++) {
		cout << rng.inRange(10, 5) <<"\t";
	}

	return 0;
}

BasicLCG::BasicLCG()
{
	iCurrent = 0; //Chose a default seed 
}

BasicLCG::BasicLCG(unsigned long iSeed)
{
	iCurrent = iSeed;
}

void BasicLCG::seed(unsigned long iSeed)
{
	iCurrent = iSeed;
}

//NOTE: a and c are selected for no particular properties
// and I have not run statistical tests on this number gernerator
// it is was written for demonstration purposes only.
unsigned long BasicLCG::nextNumber()
{
	unsigned long iOutput;
	unsigned long iTemp;
	int i;
	//I only want to take the top two bits
	//This will shorten our period to (2^32)/16=268,435,456
	//Which seems like plenty to me.
	for(i=0; i<16; i++)
	{
		//Since this is mod 2^32 and our data type is 32 bits long
		// there is no need for the MOD operator.
		iCurrent = (3039177861 * iCurrent + 1);
		iTemp = iCurrent >> 30;
		iOutput = iOutput << 2;
		iOutput = iOutput + iTemp;
	}
	return iOutput;	
}

unsigned short int BasicLCG::nextInt()
{
	unsigned short int iOutput;
	unsigned long iTemp;
	int i;
	//No need to limit ourselves...
	for(i=0; i<8; i++)
	{
		//Since this is mod 2^32 and our data type is 32 bits long
		// there is no need for the MOD operator.
		iCurrent = (3039177861 * iCurrent + 1);
		iTemp = iCurrent >> 30;
		iOutput = iOutput << 2;
		iOutput = iOutput + (short int)iTemp;
	}
	return iOutput;	
}

unsigned char BasicLCG::nextChar()
{
	unsigned char cOutput;
	unsigned long iTemp;
	int i;
	for(i=0; i<4; i++)
	{
		iCurrent = (3039177861 * iCurrent + 1);
		iTemp = iCurrent >> 30;
		cOutput = cOutput << 2;
		cOutput = cOutput + (char)iTemp;
	}
	return cOutput;	
}

int BasicLCG::nextBit()
{
	iCurrent = (3039177861 * iCurrent + 1);
	return iCurrent >> 31;	 
}

double BasicLCG::nextDouble()
{
	return (double)nextNumber()/0xFFFFFFFF;
}

int BasicLCG::inRange(int iMin, int iMax)
{
	int Diff;
	//IF the user put them in backwards then swap them
	if (iMax<iMin)
	{ 
		//Integer swap
		iMax = iMax ^ iMin; //iMax holds iMax ^ iMin
		iMin = iMax ^ iMin; //iMin= (iMax ^ iMin) ^ iMin = iMax (original)
		iMax = iMax ^ iMin; //iMax= (iMax ^ iMin) ^ iMax = iMin (original)
	}
	Diff = iMax - iMin + 1;
	return (int) (nextDouble()*Diff)+iMin;
}


IV Topics
*This is just a side note. I know that slot machines and online games tend use LCG's to do thier random numbers and I wonder if the LCGs' inherent patterns arn't senced by the people who play these games. In MMRPG's people have 'superstitions' about how to do this or that based upon patterns... very often these people seem to have a statistical advantage (though that is not a scientific statment). Could it be that these users have found a way to catch a bit wave? to sync into a pattern in the LCG? My brain tells me that those random number generators are working pretty hard to keep up with the demand of the users and 1 user falling into a pattern seems unlikely.

V Referances
Knuth, Donald E., The Art of Computer Programming V2 Seminumerical Algorithms 3rd Ed., Addison-Wesley, pp10-11, 170-173, 1997.

Stephen K. Park and Keith W. Miller, Random Number Generators: Good Ones Are Hard To Find, Communications of the ACM, 31(10):1192-1201, 1988.

Linear congruential generator. Wikipedia, http://en.wikipedia....ntial_generator


-------------------------------------------------------------------
Hope someone finds this interesting/helpful

If there is interest next time I will tackle the Mersenne twister which is a faster slightly more secure algorithm.

This post has been edited by NickDMax: 16 June 2009 - 08:29 PM


Is This A Good Question/Topic? 1
  • +

Replies To: Random Number Generation 102

#2 k0b13r  Icon User is offline

  • D.I.C Head
  • member icon

Reputation: 15
  • View blog
  • Posts: 243
  • Joined: 18-July 06

Posted 16 April 2007 - 12:35 AM

Very nice tutorial :D
Was This Post Helpful? 0
  • +
  • -

#3 Suhair  Icon User is offline

  • New D.I.C Head

Reputation: 0
  • View blog
  • Posts: 1
  • Joined: 16-June 07

Posted 16 June 2007 - 09:50 PM

Hello,
Thank you for your code and I have another question. How can we write a program to get a , b if I have xi

Bye
Was This Post Helpful? 0
  • +
  • -

#4 NickDMax  Icon User is offline

  • Can grep dead trees!
  • member icon

Reputation: 2250
  • View blog
  • Posts: 9,245
  • Joined: 18-February 07

Posted 21 July 2007 - 09:47 PM

Quote

How can we write a program to get a , b if I have xi?


Lets see... what are the a, and c and M needed to generate 133456? Who knows!!!

Ok, but what if I have a sequence of numbers Xi how can I figure out the LCG used?

If you play with LCG's enough you begin to realize that it is not hard to statistically guess at the next number.

Now at the moment I don't have any reference material to look over, or a compiler to work with but here are my thoughts:

There are only a few LCG's commonly used out there, so if you have access to the random number generator and can seed it with a known value it really is not hard to determine which LCG is used. OF course it would be wishful thinking to imagine that we have such access. If you have a large sample of a generated sequence then you can generate a statistical profile that will narrow the list to just a few which can then be checked.

I don't know if it is algebraically possible to solve for the LCG used, but I imagine that it is possible to quickly narrow down the possibilities. The LCG's have restrictions which should allow you to look for an intersection of sets as the size of the sequence expands.

I know that it is statistically possible to guess the next number for a given LCG, and this basically means that you are guessing at LCG's and should be able to narrow them down to one.

I think there is an exercise in Knuth's book on this vary matter.
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1