#include "common.h"
#include "Halton.h"



void Halton::TestHalton(int n, int dim)
{
	Halton h(dim);

	for (int i = 1; i <= n; ++ i)
	{
		std::cout << "halton " << i << " of dim " << dim << " = " << h.GetNext() << std::endl;
	}
}


void Halton::TestPrime()
{
	int dim = 2;

	for (int i = 1; i < 10; ++ i)
	{
		std::cout << "prime for dim " << i << " " << Halton::FindPrime(i) << std::endl;
	}
}


Halton::Halton(int dim): mIndex(0)
{
	mBase = FindPrime(dim);
}


float Halton::GetNext()
{
	++ mIndex;

	float result = .0f;
	float fraction = 1.0f / (float)mBase;
	int idx = mIndex;

	while (idx > 0) 
	{
		int digit = idx % mBase;
		result += fraction * (float)digit;

		idx  = (idx - digit) / mBase;
		fraction /= (float)mBase;
	}

	//mIndex %= 100000000;
	return result;
}


bool Halton::IsPrime(int n) 
{
	bool isPrime = true;

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

	return isPrime;
}
 

int Halton::FindPrime(int idx) 
{
	if (idx < 1) 
	{
		std::cerr << "error: cannot find " << idx << "th prime" << std::endl;
		return 0;
	}

	// only even prime numbers
	if (idx == 1) return 2;

	int number = 3;
	int primeFound = 1;

	while (1)
	{
		if (IsPrime(number)) ++ primeFound;
		if (primeFound == idx) break;

		number += 2;
	}

	return number;
}


void HaltonSequence::TestHalton(int n, int dim)
{
	HaltonSequence h(dim);
	float *haltons = new float[dim];

	for (int i = 1; i <= n; ++ i)
	{
		h.GetNext(haltons);
		std::cout << "halton " << i << " of dim " << dim << " = ";

		for (int j = 0; j < dim; ++ j)
			std::cout << haltons[j] << " ";

		std::cout << std::endl;
	}
}


HaltonSequence::HaltonSequence(int dim)
{
	for (int i = 0; i < dim; ++ i)
	{
		mHaltons.push_back(Halton(i + 1));
	}
}


void HaltonSequence::GetNext(float *numbers)
{
	for (size_t i = 0; i < mHaltons.size(); ++ i)
	{
		numbers[i] = mHaltons[i].GetNext();
	}
}