#ifndef __HALTON_H
#define __HALTON_H

#include <iostream>
using namespace std;

namespace GtpVisibilityPreprocessor {

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

template<int T>
class Halton {
  float _invBases[T];
  float _prev[T];
  
public:
  
  void Reset() {
	for (int i=0; i < T; i++) 
	  _prev[i] = 0;
  }
  
  Halton() {
	for (int i=0; i < T; i++) {
	  int base = FindPrime(i+1);
	  if (base == 1)
		base++;
	  _invBases[i] = 1.0f/base;
	}
	Reset();
  }
  
  void
  GetNext(float *a) {
	for (int i=0; i < T; i++) {
	  a[i] = halton(_invBases[i], _prev[i]);
	  _prev[i] = a[i];
	}
  }
  
};

class 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;
  }
};


/**
 * 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) {
  //  if (index < 1) {
  //	cerr<<"FindPrime: The argument must be non-negative."<<endl;
  //	return -1;
  //  }

  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;
}

struct HaltonSequence {
public:
  int index;

  HaltonSequence():index(1) {}

  void Reset() {
	index = 1;
  }

  void GenerateNext() {
	index++;
  }
  
  /**
   * 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 GetNumber(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 {
	  cerr<<"Error generating Halton sequence."<<endl;
	  exit(1);
	}
  }
};

extern Halton2 halton2;
}

#endif