mlib - multimedia libraries Code

#include <iostream>
#include <cmath>
#include <time.h>
#include "random.h"
using namespace std;

#pragma warning(disable: 4222)

int random::_n = 0;

random::random()
{
	init();
	set_seed();
}

random::random(long seed)
{
	init();
	set_seed(seed);
}

random::~random()
{
	cleanup();
}

// result in [begin, end)
int random::integer(const int begin, const int end)
{
	return int(floor(begin + (end - begin)*ran2(&_seed)));
	/*
	float r = ran2(&_seed);
	cout << "r=" << r << endl;
	double z = begin + (end - begin)*r;
	cout << "z=" << z << endl;
	return int(floor(z));
	*/
}

double random::real(const double begin, const double end)
{
	return begin + (end - begin)*double(ran2(&_seed));
}

/*
NR function gasdev()
Returns a normally distributed deviate with zero mean and unit variance, using ran2()
as the source of uniform deviates.
*/
double random::normal(double mean, double sd)
{
	static bool iset = false;
	static float gset;
	double rsq, v1, v2;

	if(!iset) {		// We don't have an extra deviate handy, so
		do {
			v1 = 2.0*ran2(&_seed) - 1.0;	// pick two uniform numbers in the square extending
			v2 = 2.0*ran2(&_seed) - 1.0;	// from -1 to +1 in each direction,
			rsq = v1*v1 + v2*v2;			// see if they are in the unit circle,
		} while(rsq >= 1.0 || rsq == 0.0);				// and if they are not, try again.
		const float fac = sqrt(-2.0*log(rsq)/rsq);
		// Now make the Box-Muller transformation to get two normal deviates. Return one and save the other for next time.
		gset = v1*fac;
		iset = true;			// Set flag.
		return v2*fac*sd + mean;
	}
	else {						// We have an extra deviate handy,
		iset = false;			// so unset the flag,
		return gset*sd + mean;	// and return it.
	}
}

void random::shuffle(void *p, int n, unsigned int stride)
{
	char *temp = new char[stride];
	char *zero = (char *) p;
	for(int i = 0, j; i < n; i++) {
		j = integer(i, n);
		memcpy(temp, zero + stride*i, stride);
		memcpy(zero + stride*i, zero + stride*j, stride);
		memcpy(zero + stride*j, temp, stride);
	}
	delete [] temp;
}

void random::permutation(int n, vector<int> &out)
{
	out.clear();
	for(int i = 0; i < n; i++) out.push_back(i);
	shuffle(out);
}

void random::set_seed()
{
	set_seed(time(NULL) + 100*_n++);
}

void random::set_seed(long seed)
{
	if(seed > 0)
		_seed = -seed;
	else
		_seed = seed;
}

long random::get_seed()
{
	return _seed;
}

#define RANDOM_IM1		2147483563
#define RANDOM_IM2		2147483399
#define RANDOM_AM		(1.0f/RANDOM_IM1)
#define RANDOM_IMM1		(RANDOM_IM1-1)
#define RANDOM_IA1		40014
#define RANDOM_IA2		40692
#define RANDOM_IQ1		53668
#define RANDOM_IQ2		52774
#define RANDOM_IR1		12211
#define RANDOM_IR2		3791
#define RANDOM_NTAB		32
#define RANDOM_NDIV		(1 + RANDOM_IMM1/RANDOM_NTAB)
#define RANDOM_EPS		1.2e-7
#define RANDOM_RNMX		(1.0f - float(RANDOM_EPS))

/*
Long period (> 2 × 1018) random number generator of L’Ecuyer with Bays-Durham shuffle
and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of
the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter
idum between successive deviates in a sequence. RANDOM_RNMX should approximate the largest floating
value that is less than 1.
*/
float random::ran2(long *idum)
{
	int j;
	long k;
	float temp;

	if(*idum <= 0) {		// Initialize.
		if(-(*idum) < 1)
			*idum = 1;		// Be sure to prevent idum = 0.
		else
			*idum = -(*idum);
		idum2 = (*idum);
		for(j = RANDOM_NTAB + 7; j >= 0; j--) {	// Load the shuffle table (after 8 warm-ups).
			k = (*idum)/RANDOM_IQ1;
			*idum = RANDOM_IA1*(*idum - k*RANDOM_IQ1) - k*RANDOM_IR1;
			if(*idum < 0)
				*idum += RANDOM_IM1;
			if(j < RANDOM_NTAB)
				iv[j] = *idum;
		}
		iy = iv[0];
	}
	k = (*idum)/RANDOM_IQ1;						// Start here when not initializing.
	*idum = RANDOM_IA1*(*idum - k*RANDOM_IQ1) - k*RANDOM_IR1;	// Compute idum=(RANDOM_IA1*idum) % RANDOM_IM1 without overflows by Schrage’s method.
	if(*idum < 0)
		*idum += RANDOM_IM1;
	k = idum2/RANDOM_IQ2;
	idum2 = RANDOM_IA2*(idum2 - k*RANDOM_IQ2) - k*RANDOM_IR2;	// Compute idum2=(RANDOM_IA2*idum) % RANDOM_IM2 likewise.
	if(idum2 < 0)
		idum2 += RANDOM_IM2;
	j = iy/RANDOM_NDIV;							// Will be in the range 0..RANDOM_NTAB-1.
	iy = iv[j] - idum2;							// Here idum is shuffled, idum and idum2 are combined to generate output.
	iv[j] = *idum;
	if(iy < 1)
		iy += RANDOM_IMM1;
	if((temp = RANDOM_AM*iy) > RANDOM_RNMX)
		return RANDOM_RNMX;						// Because users don’t expect endpoint values.
	else
		return temp;
}

void random::init()
{
	idum2 = 123456789;
	iy = 0;
	iv = new long[RANDOM_NTAB];
}

void random::cleanup()
{
	delete[] iv;
	iv = 0;
}