Slumptalsgenerator (WELL)

/**
 * \class	Random
 * 
 * \ingroup	core
 * \brief	Random number generator
 * \author	Mattias Gustavsson	
 * 
 * Random number generation, using the WELL algorithm by F. Panneton, P. L'Ecuyer and M. Matsumoto.
 * The creators of the algorithm describes it like this:
 *
 *      Fast uniform random number generators with extremely long periods have been defined and implemented based on 
 *      linear recurrences modulo 2. The twisted GFSR and the Mersenne twister are famous recent examples. Besides the 
 *      period length, the statistical quality of these generators is usually assessed via their equidistribution 
 *      properties. The huge-period generators proposed so far are not quite optimal in that respect. In this paper, 
 *      we propose new generators, with better equidistribution and "bit-mixing" properties for equivalent period length 
 *      and speed. Approximately half of the coefficients of the characteristic polynomial of these generators are 
 *      nonzero. The state of our new generators evolves in a more chaotic way than for the Mersenne twister. We 
 *      illustrate how this can reduce the impact of persistent dependencies among successive output values, which can 
 *      be observed in certain parts of the period of gigantic generators such as the Mersenne twister.
 *
 * More information in the original paper: http://www.iro.umontreal.ca/~panneton/WELLRNG.html
 *
 * This code is originally based on WELL512 C/C++ code written by Chris Lomont (published in Game Programming Gems 7) 
 * and placed in the public domain. There's also an article about it on Lomont's site:
 *
 *      http://lomont.org/Math/Papers/2008/Lomont_PRNG_2008.pdf
 */

#ifndef Random_H
#define Random_H

class Random
	{
	public:
		/**
		 * Constructor. Initialises generator with a default seed
		 */
		Random();

		/**
		 * Seeds the generator with a custom seed value. 
		 */
		void Seed(
			unsigned int seed	///< Custom seed value
			);

		/**
		 * Retrieves the internal state of the generator. This consists of 16
		 * unsigned ints, and one index value. This state can be stored and
		 * later re-applied to a generator to continue along the same sequence.
		 */
                void GetState(
			unsigned int state[16], ///< Array to receive the internal state from the generator
			int& index ///< Variable to receive the index value from the generator
			);

		/**
		 * Re-applies a state that was previously retrieved by calling GetState.
		 * Useful to continue the generator from a previously saved point.
		 */
                void SetState(
			const unsigned int state[16], ///< Array of state values to be applied to the generator
			int index // Index value to be applied to the generator
			);

		/** 
		 * Generates a random number on [0,0xffffffff]-interval
		 */
		unsigned int RandomInteger();

		/**
		 * Generates a random number on [0,0.99999999999...] interval
		 */
		float RandomFloat();

		/**
		 * Generates a random number on [min,max] interval
		 */
		int RandomRanged(int min, int max);

	private:
                unsigned int state_[16];	///< The internal state of the generator
                unsigned int index_; ///< Current index value

};

#endif /* Random_H */

//*** Random.cpp ***

#include "Random.h"

//*** Constructor ***

Random::Random():
        index_(0)
	{
        Seed(0x5ee39c34);  // Default seed
        }


//*** Seed ***

void Random::Seed(unsigned int s)
	{	
	// Some bit manipulation magic to propagate seed value to the internal state
        state_[0]=s^0xf68a9fc1;
        for (int i=1; i<16; i++) 
		{
                state_[i] = (0x6c078965U * (state_[i-1] ^ (state_[i-1] >> 30)) + i); 
		}
	index_ = 0;
	}


//*** SetState ***

void Random::SetState(const unsigned int state[16], int index)
        {
        for (int i=0; i<16; i++)
            {
            state_[i]=state[i];
            }

        index_=index&15;
        }


//*** GetState ***

void Random::GetState(unsigned int state[16], int& index)
        {
        for (int i=0; i<16; i++)
            {
            state[i]=state_[i];
            }

        index=index_;
        }


//*** RandomInteger ***

unsigned int Random::RandomInteger()
	{
	// This is basically just the WELL algorithm as-is
        unsigned int a = state_[index_];
        unsigned int c = state_[(index_+13)&15];
        unsigned int b = a^c^(a<<16)^(c<<15);
        c = state_[(index_+9)&15];
        c ^= (c>>11);
        a = state_[index_] = b^c;
        unsigned int d = a^((a<<5)&0xda442d24U);
        index_ = (index_ + 15)&15;
        a = state_[index_];
        state_[index_] = a^b^d^(a<<2)^(b<<18)^(c<<28);
        return state_[index_];
	}


//*** RandomRanged ***

int Random::RandomRanged(int min, int max)
	{
	// Calculate the range. The +1 is because the RandomFloat method returns
	// values up to, but not including, 1. So, a max/min difference of 0, will 
	// result in a range value of 1, which will (correctly) return all zeros.
	// If we didn't do this, probability that max would come up would be wrong.
	int range=(max-min)+1;

	// If min was greater than max, just return min
	if (range<=0)
		{
		return min;
		}
	
	// Get a random float and scale it with the range. By generating a float and
	// scaling it by range, we guarantee an equal distribution of values, which
	// we are not sure to get if using modulo operator for the range.
	int value=(int)(RandomFloat()*range);
        return min+value; 
	}


//*** RandomFloat ***

float Random::RandomFloat()
	{
        // Get a random integer, and divide by 2^32
        return (RandomInteger()/4294967296.0f);     
	}