#include "Random.h"
#include <limits>
#include <iostream> // debug

extern unsigned long long cntPow; // debug
//extern std::ofstream dbg; // debug

#include "zetaArray.cpp"

double Urnd::rnd() {
	static std::random_device rd;
	static std::mt19937_64 gen(rd()); // Mersenne twister (64 bit)
	static std::uniform_real_distribution<double> urd(0., 1.);
	return urd(gen);
}
int Irnd::rnd(int low, int high) {
	static std::random_device rd;
	static std::mt19937_64 gen(rd()); // Mersenne twister (64 bit)
	using Dist = std::uniform_int_distribution<int>;
	static Dist uid{};
	return uid(gen, Dist::param_type{ low, high });
}
int Grnd::rnd(double p) {
	Urnd urnd;
	double U = urnd.rnd();
	int G = static_cast<int>(floor(log(1 - U) / log(p)));
	return G;
}
int Prnd::rndslow() {
	double const beta = 1.5;
	double const zeta = // zeta(3/2)
		2.61237534868548834334856756792407163057080065240006340757332824881492776768827286099624386812631195238297;
	Urnd urand;
	double U = 2 * urand.rnd() - 1;
	double target = std::abs(U) * zeta;
	// slow method
	double sum = 1.;
	unsigned long k = 1;
	while (target > sum) {
		k += 1;
		sum += pow(static_cast<double>(k), -beta);
		cntPow++;
		if (k >= 30000) break; // ### FAKE -> must be designed and programmed properly
	}
	
	if (k > INT32_MAX) k = INT32_MAX;
	int P = static_cast<int>(k);
	if (U < 0) P = -P;
	return P;
}

inline int sgn(double x) {
	return (x < 0) ? -1 : 1;
}

// Function to find insert position of K in sorted array
static int findIndex(double arr[], int n, double K)
{
	// Lower and upper bounds
	int start = 0;
	int end = n - 1;
	// Traverse the search space
	while (start <= end) {
		int mid = (start + end) / 2;
		// If K is found
		if (arr[mid] == K)
			return mid;
		else if (arr[mid] < K)
			start = mid + 1;
		else
			end = mid - 1;
	}
	// Return insert position
	return end + 1;
}

int Prnd::rnd() {
	double const beta = 1.5;
	double const zeta = // zeta(3/2)
		2.61237534868548834334856756792407163057080065240006340757332824881492776768827286099624386812631195238297;
	unsigned int const maxIndex = 99999;

	Urnd urand;
	double U = 2 * urand.rnd() - 1;
	double target = std::abs(U) * zeta;
	if (target <= zetaArray[maxIndex]) {
		// find smallest index with zetaArray[index] > target
		int index = findIndex(zetaArray, maxIndex + 1, target);
		return (index + 1) * sgn(U);
		//int a = 0, b = maxIndex;
		//int loopcnt = 0;
		//do {
		//	if (zetaArray[a] > target) return (a + 1) * sgn(U);
		//	int c = (a + b) / 2;
		//	if (target == zetaArray[c]) return (c + 1) * sgn(U);
		//	if (target > zetaArray[c]) {
		//		a = c + 1;
		//	}
		//	else {
		//		b = c;
		//	}
		//	loopcnt++;
		//	if (loopcnt > 20) {
		//		std::cout << target << " " << a << " " << b << " " << c << std::endl;
		//	}
		//} while (a <= b);

		//throw;
		//return -1;
	}
	else { // target > zetaArray[maxIndex]

		if (target < zetaArray1[0]) { // k is between 100000 and 200000
			double sum = zetaArray[maxIndex]; // k starts with 100000
			int k = 100000;
			while (target > sum) {
				k += 1;
				sum += pow(static_cast<double>(k), -beta);
			}
			return k * sgn(U);
		}
		else {
			// find smallest index with zetaArray1[index] > target
			int index = findIndex(zetaArray1, 21400, target);

			//int a = 0, b = 21399, ti;
			//do {
			//	if (zetaArray1[a] > target) { ti = (a + 1); break; }
			//	int c = (a + b) / 2;
			//	if (target == zetaArray1[c]) { ti = (c + 1); break; }
			//	if (target > zetaArray1[c]) {
			//		a = c + 1;
			//	}
			//	else {
			//		b = c;
			//	}
			//} while (a <= b);

			double sum = zetaArray1[index];
			int k = index * 100000+100000;
			while (target > sum) {
				k += 1;
				sum += pow(static_cast<double>(k), -beta);
			}
			return k * sgn(U);
		}

		return INT32_MAX * sgn(U);
	}
}

#include <iostream>
#include <fstream>
#include <iomanip>
void Prnd::makeArray() {
	double const beta = 1.5;
	double const zeta = // zeta(3/2)
		2.61237534868548834334856756792407163057080065240006340757332824881492776768827286099624386812631195238297;
	double sum = 1.;	
	std::ofstream out;
	out.open("array.txt");
	out << "1.,";

	for (int i = 2; i <= 100000; i++) {
		sum += pow(static_cast<double>(i), -beta);
		out << std::setprecision(10) << sum << ",";
		if (i % 10 == 0) out << std::endl;
	} 
	out.close();
	std::cout << sum / zeta << std::endl;
	
}
void Prnd::makeArray2() {
	double const beta = 1.5;
	double const zeta = // zeta(3/2)
		2.61237534868548834334856756792407163057080065240006340757332824881492776768827286099624386812631195238297;
	double sum = zetaArray[99999]; // starting point at index = 100000
	std::ofstream out;
	out.open("array1.txt");

	for (int i = 1; i <= 21400; i++) {
		for (int j = 1; j <= 100000; j++) {
			sum += pow(static_cast<double>(i * 100000 + j), -beta);
		}
		out << std::setprecision(10) << sum << ",";
		if (i % 10 == 0) out << std::endl;		
	}
	out.close();
	std::cout << sum / zeta << std::endl;

}
