#include <iostream>
#include <omp.h>
#include <cstdlib>
#include <ctime>
#include <vector>

const double dx = 0.0000001;
const int monteCarloCount = 900000;

double function(double x) {
	return 4.0 / (1 + x*x);
}

void SingleThread() {
	long lLength = 1 / dx;
	double dSum = 0;
	for (long x = 0; x<lLength; ++x)
		dSum += function(x * dx) * dx;

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

void ThreadedWithCommon() {
	long lLength = 1 / dx;
	double adSum[8] = { 0,0,0,0,0,0,0,0 };
#pragma omp parallel
	{
		const int nThreadNum = omp_get_thread_num();
		const int nTotalThreads = omp_get_num_threads();
		for (long x = nThreadNum; x<lLength; x += nTotalThreads)
			adSum[nThreadNum] += function(x * dx) * dx;
	}
	double dSum = 0;
	for (int i = 0; i < 8; ++i) {
		dSum += adSum[i];
	}

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

void Threaded() {
	long lLength = 1 / dx;
	double adSum[8] = { 0,0,0,0,0,0,0,0 };
#pragma omp parallel
	{
		double dSum = 0;
		const int nThreadNum = omp_get_thread_num();
		const int nTotalThreads = omp_get_num_threads();
		for (long x = nThreadNum; x<lLength; x += nTotalThreads)
			dSum += function(x * dx) * dx;
		adSum[nThreadNum] = dSum;
	}
	double sum = 0;
	for (int i = 0; i < 8; ++i) {
		sum += adSum[i];
	}
	std::cout << sum << std::endl;
}

void ThreadedForReduction() {
	long lLength = 1 / dx;
	double dSum = 0;
#pragma omp parallel for reduction (+:dSum)
	for (long x = 0; x<lLength; ++x)
		dSum += function(x * dx) * dx;
	std::cout << dSum << std::endl;
}

void MonteCarloSingleThread() {
	int in = 0;
	double x, y;
	for (int i = 0; i < monteCarloCount; ++i) {
		x = rand() / (double)RAND_MAX;
		y = rand() / (double)RAND_MAX;
		if (x*x + y*y<1)
			++in;
	}
	std::cout << 4.0*in / monteCarloCount << std::endl;
}

void MonteCarloThreaded() {
	const int n = 1;
	std::cout << n << " threads" << std::endl;
	std::vector<int> in_total;
	unsigned int in_sum = 0;
	//#pragma omp parallel num_threads(n)
	{
		srand(int(time(NULL)) ^ omp_get_thread_num());
		unsigned int in = 0, iterations = monteCarloCount / n;
		double x, y;
		for (int i = 0; i < iterations; ++i) {
			x = rand() / (double)RAND_MAX;
			y = rand() / (double)RAND_MAX;
			if (x*x + y*y < 1)
				++in;
		}
		in_total.push_back(in);
	}
	for (int i = 0; i < in_total.size(); i++) {
		in_sum += in_total.at(i);
	}
	std::cout << 4.0*in_sum / monteCarloCount << std::endl;
}

struct StopWatch {
	int nClockStart;
	int nClockEnd;
	void Start() {
		nClockStart = clock();
	}
	void End() {
		nClockEnd = clock();

		double dDiffTicks = nClockEnd - nClockStart;
		double dDiffTime = (dDiffTicks) / (CLOCKS_PER_SEC / 1000);
		std::cout << "Time: " << dDiffTime << "ms" << std::endl;
	}
};

int main(int argc, char *argv[]) {
		StopWatch stopWatch;

		std::cout << "Single Thread" << std::endl;
		stopWatch.Start();
		SingleThread();
		stopWatch.End();


		std::cout << "Threaded with common variable" << std::endl;
		stopWatch.Start();
		ThreadedWithCommon();
		stopWatch.End();

		std::cout << "Threaded with variable in each thread" << std::endl;
		stopWatch.Start();
		Threaded();
		stopWatch.End();

		std::cout << "Threaded with for reduction" << std::endl;
		stopWatch.Start();
		ThreadedForReduction();
		stopWatch.End();

		std::cout << "Monte Carlo Threaded: ";
		stopWatch.Start();
		MonteCarloThreaded();
		stopWatch.End();

		std::cout << "Monte Carlo Single Thread" << std::endl;
		stopWatch.Start();
		MonteCarloSingleThread();
		stopWatch.End();

		std::cin.get();

	return 0;
}