random.h

随机数生成源代码

/*

	random.h
	copyright (c) 2002
	Copious Systems

*/



#ifndef _H_RANDOM_COPIOUS_2002_
#define _H_RANDOM_COPIOUS_2002_



#include <string>
using namespace std;

// In general a generator is an object that emits numbers via
// a value function.

class I_generator {
public:
	virtual void init_generator(unsigned long int *s) {};
};

class generator : public I_generator {
public:

	generator() {
					i0 = 0;
					s0 = &i0;
					init_generator(s0);
				}

	generator(unsigned long int i) {
					i0 = i;
					s0 = &i0;
					init_generator(s0);
				}

	virtual ~generator() {}

protected:

	unsigned long int *s0;
	unsigned long int i0;

//  Add a public method value() returning template determined type
	
public:

	inline unsigned long int i_stored() { return(i0); }


};


//
//	psuedo_random
//
//	From Schrage's algorithm for approximate factorization of 
//	a prime modulus,  m,
//		
//		m = a*q + r
//
//  Then a psuedo generation squence, e.g. I(j+1) = I(j) mod m,
//	can be generated without over flow. Where z in (0,m-1),
//	it can be shown that a(z mod q) and (r[z/q]) are both in (0,m-1).
//	(See Press, et. al pg. 282).
//	
//
//	We can constuct a psuedo random generator object for use in 
//  algorithms creating very long sequences.
//
//
//



template<class T,class C>
class psuedo_random : public generator, public C {

public:

	psuedo_random() {
	}

	psuedo_random(unsigned long int i) : generator(i) {
	}

	virtual ~psuedo_random() {}

	virtual void init_generator(unsigned long int *s)  {
		// A generator with no initialization.
	}

	void reinit(unsigned long int i) { i0 = i; }


	inline T m_value() { // Using the mask
		T result;
		unsigned long int k;
		long int i;
	//-----
		i0 ^= MASK();
		k = i0/Q();

		i = A()*(i0 - k*Q()) - R()*k; // Shrage's method.
		if ( i < 0 ) i += M();

		result = SAMPLE()*i; // Contract by sample spacing.

		i0 = i;
		i0 ^= MASK();
	//-----
		return(result);
	}

	inline T value() { // not using the mask.
		T result;
		unsigned long int k;
		long int i;
	//-----
		k = i0/Q();

		i = A()*(i0 - k*Q()) - R()*k; // Shrage's method.
		if ( i < 0 ) i += M();

		result = SAMPLE()*i; // Contract by sample spacing.

		i0 = i;
	//-----
		return(result);
	}


	inline unsigned long int i_value() {
		unsigned long int k;
		long int i;
	//-----
		k = i0/Q();

		i = A()*(i0 - k*Q()) - R()*k; // Shrage's method.
		if ( i < 0 ) i += M();

		i0 = i;

		return(i);
	}

};



template<class T>
class constants_1 {
public:
	constants_1() {}

	inline const unsigned long int M() { return(2147483647); }
	inline const unsigned long int A() { return(16807); }
	inline const unsigned long int Q() { return(127773); }
	inline const unsigned long int R() { return(2836); }
	inline const T SAMPLE()	{
								const T SAMPLE = T( 1.0 / T(M()) );
								return(SAMPLE);
							};
	inline const unsigned long int MASK() { return(123459876); }
};

template<class T>
class constants_r2_1 {
public:
	constants_r2_1() {}

	inline const unsigned long int M() { return(2147483563); }
	inline const unsigned long int A() { return(40014); }
	inline const unsigned long int Q() { return(53668); }
	inline const unsigned long int R() { return(12211); }
	inline const T SAMPLE()	{
								const T SAMPLE = T( 1.0 / T(M()) );
								return(SAMPLE);
							};
	inline const unsigned long int MASK() { return(123459876); }
};

template<class T>
class constants_r2_2 {
public:
	constants_r2_2() {}

	inline const unsigned long int M() { return(214748339); }
	inline const unsigned long int A() { return(40692); }
	inline const unsigned long int Q() { return(52774); }
	inline const unsigned long int R() { return(3791); }
	inline const T SAMPLE()	{
								const T SAMPLE = T( 1.0 / T(M()) );
								return(SAMPLE);
							};
	inline const unsigned long int MASK() { return(123459876); }
};


#define psuedoR1(T) psuedo_random<T, constants_1<T> >

#define psuedoR21(T) psuedo_random<T, constants_r2_1<T> >
#define psuedoR22(T) psuedo_random<T, constants_r2_2<T> >


template<class T>
class shuffle_constants {
public:
	shuffle_constants() {}
	inline const unsigned int TABLE_SIZE() { return(32); }
	inline const unsigned long int LOG2_TABLE_SIZE() { return(5); }

	inline const T EPS()	{
									const T EPS = 3.0e-16;
									return(EPS);
								}

	inline const T RNMX()	{
									const T RNMX = 1.0 - EPS();
									return(RNMX);
								}

};



template<class T>
class generator_parameters {
public:

	generator_parameters() : y(0), buffer(NULL) {
		lb = T(0); ub = T(0); span = T(0);
	}

	virtual ~generator_parameters() {}

protected:

	unsigned long int y;
	unsigned long int *buffer;

	T lb, ub, span;

public:

	void setBounds(T lower, T upper) {
		if ( lower < upper ) {
			lb = lower;
			ub = upper;
		} else {
			ub = lower;
			lb = upper;
		}
		span = ub - lb;
	}


	inline T shift(T &result) {
		result *= span;
		result += lb;
		return(result);
	}

};


//
//
// random_generator
//


template<class T, class C>
class random_generator : public generator, public C, public generator_parameters<T>  {
public:

	random_generator() : generator_parameters<T>(), psr(0), generator() {
		random_generator::init_generator(s0);
	}

	random_generator(int i) : generator_parameters<T>(), psr(i), generator(1) {
		random_generator::init_generator(s0);
	}

	virtual ~random_generator() { }

	inline const unsigned long int NDIV() {
		const unsigned long int MSHIFT = 1 + (psr.M()>>LOG2_TABLE_SIZE());
		return(MSHIFT);
	}
	inline const unsigned long int TABLE_CONTRACTION() { return(y/NDIV()); }


private:

	psuedoR1(T) psr;

	virtual void init_generator(unsigned long int *s)
		{
			if ( buffer == NULL ) buffer = new unsigned long int (TABLE_SIZE()+10);
			if ( buffer == NULL ) throw "Bad Allocation";

			if ( *s == 0 ) *s = 1;
			for ( unsigned int j = TABLE_SIZE() + 7; j > 0; j-- ) {
				unsigned long int i = psr.i_value();
				if ( j < TABLE_SIZE() ) buffer[j] = i;
			}
			y = buffer[0];
		}


public:

	void reinit(unsigned long int i) { i0 = i; init_generator(s0); }

	inline T value() { // not using the mask.

		T result = psr.value();

		unsigned long int i = psr.i_stored();
		unsigned long int j = TABLE_CONTRACTION();

		y = buffer[j];
		buffer[j] = i;

		if ( result > RNMX() ) result = RNMX();
		return(result);
	}


	inline T pick_a_number_between_bounds() {
		T result = value();
		return(shift(result));
	}

};



#define rgen(T) random_generator<T, shuffle_constants<T> >




template<class T, class C>
class bi_random_generator : public generator, public C, public generator_parameters<T>  {
public:

	bi_random_generator() : generator_parameters<T>(), psr_1(0), psr_2(0), generator() {
		bi_random_generator::init_generator(s0);
	}

	bi_random_generator(int i) : generator_parameters<T>(), psr_1(i), psr_2(i), generator(1) {
		bi_random_generator::init_generator(s0);
	}