#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <conio.h>

#include "ComputationPI.h"
#include "mt19937.h"

// Libra
#include <Libra.h>    // C and C++ API


double arctan_pi(int termCount)
{
	gVar terms_even = div(1.f, step(1, termCount * 2 - 3, 4));
	gVar terms_odd =  div(1.f, step(3, termCount * 2 - 1, 4));

	return double(sum(terms_even - terms_odd));
}


// helper function for elementwise square -> x .* x
gVar sqr(const gVar& x)
{
	return mul(x,x);
}

double montecarlo_pi(int size, int it)
{
	int s = 0;
	for (int i = 0; i < it; i++)
	{
		gVar x = rand(size);
		gVar y = rand(size);

		s += int(sum(cast(sqr(x) + sqr(y) <= 1.f, libra_GetDefaultDataType())));
	}

	return 4.0 * s / (size * size * it);
}


int main(int argc, char** argv)
{
	std::cout << "*****************" << std::endl << std::endl;
	std::cout << "Computation of PI" << std::endl << std::endl;
	std::cout << "*****************" << std::endl << std::endl;

	if (libra_Init(argc, argv) != 0)
		return 1;

	libra_SetDefaultDataType(GFLOAT64);

	computePI();

	libra_Shutdown();
	printf("\nPress a key to exit...\n");
	_getch();
}

void computePI()
{
	// arctan(1) = PI/4.
	compute_PI_using_gregory_s_arctan_series();
	compute_PI_using_monte_carlo();
}

void compute_PI_using_gregory_s_arctan_series()
{
	std::cout << std::endl << std::endl;
	std::cout << "------------------------------------------------------" << std::endl;
	std::cout << "------- Compute PI using Gregorie's arctan series ----" << std::endl;
	std::cout << "------------------------------------------------------" << std::endl;
	std::cout << std::endl;

	const int termCount = 1000000;

	// CPU version
	{
		std::cout << std::endl << "CRT Cpu version ..." << std::endl;
		double startTime = libra_GetTime();

		double Pi = 0.0;
		for (int i=0; i < termCount; i+=2)
		{
			Pi += 1.f / (2*i + 1);
			Pi -= 1.f / (2*i + 3);
		}

		double totalTime = libra_GetTime() - startTime;

		std::cout << "Total time : " << totalTime << " seconds" << std::endl;
		std::cout << "Pi = " << Pi * 4 << std::endl;
		std::cout << std::endl;
	}

	// Libra version vectorizes code
	{
		// warmup phase to get accurate timing results.
		arctan_pi(termCount); arctan_pi(termCount); arctan_pi(termCount);

		std::cout << std::endl << "Libra version ..." << std::endl;
		double startTime = libra_GetTime();

		double Pi = arctan_pi(termCount);

		double totalTime = libra_GetTime() - startTime;
		std::cout << "Total time : " << totalTime << " seconds" << std::endl;
		std::cout << "Pi = " << Pi * 4 << std::endl;
		std::cout << std::endl;
	}
}

// mersenne twister helper function
inline float genrand_mt()
{
	return (_genrand_dc() + 1.0f) / 4294967296.0f; 
}


void compute_PI_using_monte_carlo()
{
	std::cout << std::endl << std::endl;
	std::cout << "------------------------------------------------------" << std::endl;
	std::cout << "------- Compute PI using Monte-Carlo algorithm -------" << std::endl;
	std::cout << "------------------------------------------------------" << std::endl;
	std::cout << std::endl;

	const int it = 50, size = 2048;
	const int n = it * size * size; // <-- n = number of samples

	std::cout << "using " << n << " number of samples." << std::endl;
	{
		std::cout << std::endl << "CRT Cpu version ..." << std::endl;

		double startTime = libra_GetTime();

		int s = 0;

		for (int i = 0; i < n; i++)
		{
			const float x = genrand_mt();
			const float y = genrand_mt();

			s += (x*x + y*y) <= 1.f;
		}

		double Pi = 4.0 * s / n;

		double totalTime = libra_GetTime() - startTime;

		std::cout << std::endl << "Pi = " << Pi << std::endl;
		std::cout << std::endl << "Total time : " << totalTime << " seconds" << std::endl;
		std::cout << std::endl << std::endl << std::endl;
	}
	{
		std::cout << std::endl << "Libra version ..." << std::endl;

		// warm up phase
		montecarlo_pi(size, it);

		double startTime = libra_GetTime();

		double Pi = montecarlo_pi(size, it);

		double totalTime = libra_GetTime() - startTime;

		std::cout << std::endl << "Pi = " << Pi << std::endl;
		std::cout << std::endl << "Total time : " << totalTime << " seconds" << std::endl;
		std::cout << std::endl << std::endl << std::endl;
	}
}