#ifndef RANDOM_NUMBER_GENERATORS_HPP
#define RANDOM_NUMBER_GENERATORS_HPP
#include <cmath>
#include <vector>
#include <cassert>
#include <iostream>

#define WEIGHTED_RANDOM_DEBUG 0

namespace RandomNumberGenerators {

inline float uniform()
/* Uniform random number generator x(n+1)= a*x(n) mod c
				with a = pow(7,5) and c = pow(2,31)-1.
				Copyright (c) Tao Pang 1997. */
	{
		const int ia=16807,ic=2147483647,iq=127773,ir=2836;
		int il,ih,it;
		float rc;
		static int iseed = rand();
		ih = iseed/iq;
		il = iseed%iq;
		it = ia*il-ir*ih;
		if (it > 0)
		{
			iseed = it;
		}
		else
		{
			iseed = ic+it;
		}
		rc = ic;
		return iseed/rc;
	}

inline float gaussian(float mean, float sigma)
{

	float x1, x2, w, y1, y2;
 
	do {
		x1 = 2.0 * uniform() - 1.0;
		x2 = 2.0 * uniform() - 1.0;
		w = x1 * x1 + x2 * x2;
	} while ( w >= 1.0 );

	w = sqrt( (-2.0 * log( w ) ) / w );
	y1 = x1 * w;
	y2 = x2 * w;
	
	float ret = y1*sigma + mean;

	return ret;
}

inline std::size_t uniformInteger(std::size_t upperBound=1) {

	
	/// @bug there was a man entry about how this leads to a lousy uniform
	/// @bug distribution in practice. should probably review
	assert(upperBound > 0);
	return ((rand()) % ((int)upperBound));
}

	
	

/// Randomizes from probabilistically weighted distribution. Thus,
/// sum of passed in weights should be 1.0
inline std::size_t weightedRandomNormalized(std::vector<float> weights) {

        // Choose a random bounded mass between 0 and 1
	float cutoff = ((float)(rand())) / (float)RAND_MAX;

	//std::cout << "cutoff : " << cutoff << std::endl;

    // Sum up mass, stopping when cutoff is reached. This is the typical
    // weighted sampling algorithm.
	float mass = 0;
	for (std::size_t i = 0; i< weights.size() ; i++) {
		mass += weights[i];
		//std::cout << "mass: " << mass << std::endl;
		if (mass >= cutoff)
			return i;
	}

    // Just in case something slips through the cracks
	return weights.size()-1;
}

inline std::size_t weightedRandom(const std::vector<int> & weights, unsigned int weightTotalHint = 0) {
	

	if (weightTotalHint == 0)	{		
		for (std::size_t i = 0; i < weights.size();i++) 
			weightTotalHint += weights[i];
	}
	
	const int sampledSum = uniformInteger(weightTotalHint);
	int sum = 0;
	if (WEIGHTED_RANDOM_DEBUG) std::cout << "[RNG::weightedRandom()] weightTotal = " << weightTotalHint <<
			 std::endl; 

	for (std::size_t i = 0; i < weights.size();i++) {
		if (WEIGHTED_RANDOM_DEBUG) 
		std::cout << "[RNG::weightedRandom()] weight[" << i << "] = " << weights[i] <<
			 std::endl; 

		sum += weights[i];
		if (sampledSum <= sum) {
			if (WEIGHTED_RANDOM_DEBUG) 
			std::cout << "[RNG::weightedRandom()] sampled index " << i << "(" <<
			 "running sum = " << sum << ", sampled sum = " << sampledSum << std::endl;
			return i;
		}
	}

	return weights.size()-1;
}

}
#endif