#include <iostream>
#include <cmath>
#include <random>
#include <time.h> 
using namespace std;

const double PI = 3.141592653589793;
const double PI2 = 2.0 * 3.141592653589793;
const double PI4 = 4.0 * 3.141592653589793;
const float PIf = 3.141592653589793;
const float PI2f = 2.0 * 3.141592653589793;
const float PI4f = 4.0 * 3.141592653589793;

const int N = 14400;
const int HALF_N = N >> 1;
const double STEP = 4.0 * PI / (double)N; 
const double INV_STEP = 1.0 / STEP;
static double LUT[N * 2];

const int Nsm = 90;
const int HALF_Nsm = Nsm >> 1;
const float STEPsm = 4.0 * PIf / (float)Nsm;
const float INV_STEPsm = 1.0 / STEPsm;
static float LUTsm[Nsm * 2];

int seed = (int)clock();

void benchLUT() {
	cout << "LUT (double)" << endl;
	srand(seed);
	double result = 0.0;
	clock_t t;
	t = clock();
	for(int x = 0; x < 100000; x++) {
		double r = PI4 * ((double)rand() / RAND_MAX) - PI2;
		int index = (int)((r * INV_STEP) + HALF_N) << 1;
		result += LUT[index] + LUT[index + 1];
	}
	t = clock() - t;
	cout << "Sum : " << result << endl;
	cout << "Time : " << t << " : " << ((float)t)/CLOCKS_PER_SEC << endl;
}
void benchSinCos() {
	cout << "CMATH (float)" << endl;
	srand(seed);
	float result = 0.0;
	clock_t t;
	t = clock();
	for(int x = 0; x < 100000; x++) {
		float r = (float)(PI4f * ((float)rand() / RAND_MAX) - PI2f);
		result += sin(r) + cos(r);
	}
	t = clock() - t;
	cout << "Sum : " << result << endl;
	cout << "Time : " << t << " : " << ((float)t)/CLOCKS_PER_SEC << endl;
}
void benchPoly() {
	cout << "Poly (float)" << endl;
	srand(seed);
	float result = 0.0;
	clock_t t;
	t = clock();
	for(int x = 0; x < 100000; x++) {
		float r = (float)(PI4f * ((float)rand() / RAND_MAX) - PI2f);
		// wrap
		if (r < -3.14159265) {
			r += 6.28318531;
		} else if (r > 3.14159265) {
			r -= 6.28318531;
		}
		// sin
		if (r < 0) {
			result += 1.27323954 * r + 0.405284735 * r * r;
		} else {
			result += 1.27323954 * r - 0.405284735 * r * r;
		}
		// cos
		r += 1.57079632;
		if (r > 3.14159265) {
    		r -= 6.28318531;
		}
		if (r < 0) {
    		result += 1.27323954 * r + 0.405284735 * r * r;
		} else {
    		result += 1.27323954 * r - 0.405284735 * r * r;
		}
		//result = sin(r) + cos(r);
	}
	t = clock() - t;
	cout << "Sum : " << result << endl;
	cout << "Time : " << t << " : " << ((float)t)/CLOCKS_PER_SEC << endl;
}
void benchInterp() {
	cout << "Small LUT Interp (float)" << endl;
	srand(seed);
	float result = 0.0;
	clock_t t;
	t = clock();
	for(int x = 0; x < 100000; x++) {
		float r = PI4f * ((float)rand() / RAND_MAX) - PI2f;
		int index = (int)((r * INV_STEPsm) + HALF_Nsm) << 1;
		result += LUTsm[index] + LUTsm[index + 1];
	}
	t = clock() - t;
	cout << "Sum : " << result << endl;
	cout << "Time : " << t << " : " << ((float)t)/CLOCKS_PER_SEC << endl;
}

int main() {
	// build LUT
	int i = 0;
	double r = -2.0 * PI;
	for(i=0; i < N; i++, r += STEP) {
    	LUT[(i << 1)] = sin(r);
    	LUT[(i << 1) + 1] = cos(r);
	}
	// build small LUT
	r = -2.0 * PIf;
	for(i=0; i < Nsm; i++, r += STEPsm) {
    	LUTsm[(i << 1)] = sin(r);
    	LUTsm[(i << 1) + 1] = cos(r);
	}
	
	benchLUT();
	benchSinCos();
	benchPoly();
	benchInterp();
	benchLUT();
	benchSinCos();
	benchPoly();
	benchInterp();
	return 0;
}