#ifndef __HALTON_H
#define __HALTON_H

#include <iostream>


/** Assert whether the argument is a prime number.
	@param number the number to be checked
*/
inline bool IsPrime(const int number) 
{
	bool isIt = true;

	for(int i = 2; i < number; ++ i) 
	{
		if(number % i == 0) 
		{
			isIt = false;
			break;
		}
	}

	if(number == 2) 
	{
		isIt = false;
	}

	return isIt;
}


/**
	Find the nth prime number.
	@param index the ordinal position in the sequence
*/
inline int FindPrime(const int index) 
{

	const int primes[] = {-1, 1, 3, 5, 7, 11, 13};
	if (index <= 6)
		return primes[index];

	int prime = 1;
	int found = 1;
	while(found != index) {
		prime += 2;
		if(IsPrime(prime) == true) {
			found++;
		}
	}
	return prime;
}

  
inline float halton(float baseRec, float prev) 
{
	float r = 1.0f - prev;
	
	if (baseRec < r)
		return prev + baseRec;
	
	float h = baseRec;
	
	float hh;
	
	do 
	{
		hh = h;
		h *= baseRec;
	} while (h > r);
	
	return prev + hh + h - 1.0f;
}


template<int T> struct Halton 
{
	static float _invBases[T];
	float _prev[T];

public:

	void Reset() 
	{
		for (int i=0; i < T; i++) 
			_prev[i] = 0;
	}

	Halton(const bool initializeBases) 
	{
		for (int i=0; i < T; i++) 
		{
			int base = FindPrime(i+1);
		
			if (base == 1)
				base++;
			_invBases[i] = 1.0f/base;
		}
	}

	Halton() 
	{
		Reset();
	}
  
	void GetNext(float *a) 
	{
		for (int i=0; i < T; i++) 
		{
			a[i] = halton(_invBases[i], _prev[i]);
			_prev[i] = a[i];
		}
	}

};

struct Halton2 
{
	static float _invBases[2];
	float _prev[2];

public:

	void Reset() 
	{
		_prev[0] =_prev[1] = 0;
	}

	Halton2() 
	{
		_invBases[0] = 1.0f / 2;
		_invBases[1] = 1.0f / 3;
		Reset();
	}

	void GetNext(float &a, float &b) 
	{
		a = halton(_invBases[0], _prev[0]);
		b = halton(_invBases[1], _prev[1]);

		_prev[0] = a;
		_prev[1] = b;
	}
};


  
struct HaltonSequence 
{
public:
	int index;

	static int sPregeneratedDim;
	static int sPregeneratedNumber;
	static float *sPregeneratedValues;

	// special construtor for pregenerating static halton sequences
	HaltonSequence(const int dim, const int number);

	HaltonSequence(): index(1) {}

	void Reset() 
	{
		index = 1;
	}

	void
		GetNext(const int dimensions, float *p);

	void GenerateNext() 
	{
		++ index;
	}

	double GetNumber(const int dimension)  
	{
		int base = FindPrime(dimension);
		if(base == 1) {
			base++;  //The first dimension uses base 2.
		} 

		int _p1 = base;
		float _ip1 = 1.0f/base;
		float p, u=0.0f;
		int kk, a;

		// the first coordinate
		for (p = _ip1, kk = index ; kk ;  p *= _ip1, kk /= _p1)   
			if ((a = kk % _p1))
				u += a * p;

		return u;
	}

	/**
		Returns the nth number in the sequence, taken from a specified dimension.
		@param index the ordinal position in the sequence
		@param dimension the dimension
	*/
	double GetNumberOld(const int dimension)  
	{
		int base = FindPrime(dimension);
		if(base == 1) 
		{
			++ base;  //The first dimension uses base 2.
		} 

		double remainder;
		double output = 0.0;
		double fraction = 1.0 / (double)base;
		int N1 = 0;
		int copyOfIndex = index;
		
		if ((base >= 2) && (index >= 1)) 
		{
			while(copyOfIndex > 0) 
			{
				N1 = (copyOfIndex / base);
				remainder = copyOfIndex % base;
				output += fraction * remainder;
				copyOfIndex = (int)(copyOfIndex / base);
				fraction /= (double)base;
			} 
			return output;
		}
		else 
		{
			std::cerr<<"Error generating Halton sequence."<<std::endl;
			exit(1);
		}
	}
};


#endif