#define _USE_MATH_DEFINES 1
#include <cstdio>
#include <cmath>
#include <cstdint>
#include <algorithm>
#include "windows.h"

#define PI M_PI
using int32 = int32_t;

constexpr unsigned int num_steps  = 500'000'000;
constexpr unsigned int sani_steps =   5'000'000;


constexpr double sin_garrett_c11(double x)
{
	double x2 = x * x;
	return (((((-2.05342856289746600727e-08*x2 + 2.70405218307799040084e-06)*x2
		- 1.98125763417806681909e-04)*x2 + 8.33255814755188010464e-03)*x2
		- 1.66665772196961623983e-01)*x2 + 9.99999707044156546685e-01)*x;
}

constexpr float sin_garrett_c11_s(float x)
{
	float x2 = x * x;
	return (((((-2.05342856289746600727e-08*x2 + 2.70405218307799040084e-06)*x2
		- 1.98125763417806681909e-04)*x2 + 8.33255814755188010464e-03)*x2
		- 1.66665772196961623983e-01)*x2 + 9.99999707044156546685e-01)*x;
}


constexpr float sin_spiro_s( float _fX ) {
	int32 i32I = int32( _fX * (1.0f / PI) );
	_fX = (_fX - float( i32I ) * PI);

	float fX2 = _fX * _fX;

	return (i32I & 1) ?
		-_fX * (float( 1.00000000000000000000e+00 ) +
			fX2 * (float( -1.66666671633720397949e-01 ) +
			fX2 * (float( 8.33333376795053482056e-03 ) +
			fX2 * (float( -1.98412497411482036114e-04 ) +
			fX2 * (float( 2.75565571428160183132e-06 ) +
			fX2 * (float( -2.50368472620721149724e-08 ) +
			fX2 * (float( 1.58849267073435385100e-10 ) +
			fX2 * float( -6.58925550841432672300e-13 )))))))) :
		_fX * (float( 1.00000000000000000000e+00 ) +
			fX2 * (float( -1.66666671633720397949e-01 ) +
			fX2 * (float( 8.33333376795053482056e-03 ) +
			fX2 * (float( -1.98412497411482036114e-04 ) +
			fX2 * (float( 2.75565571428160183132e-06 ) +
			fX2 * (float( -2.50368472620721149724e-08 ) +
			fX2 * (float( 1.58849267073435385100e-10 ) +
			fX2 * float( -6.58925550841432672300e-13 ))))))));
}

constexpr double sin_spiro( double _fX ) {
	int32 i32I = int32( _fX * (1.0f / PI) );
	_fX = (_fX - double( i32I ) * PI);

	double fX2 = _fX * _fX;

	return (i32I & 1) ?
		-_fX * (double( 1.00000000000000000000e+00 ) +
			fX2 * (double( -1.66666671633720397949e-01 ) +
			fX2 * (double( 8.33333376795053482056e-03 ) +
			fX2 * (double( -1.98412497411482036114e-04 ) +
			fX2 * (double( 2.75565571428160183132e-06 ) +
			fX2 * (double( -2.50368472620721149724e-08 ) +
			fX2 * (double( 1.58849267073435385100e-10 ) +
			fX2 * double( -6.58925550841432672300e-13 )))))))) :
		_fX * (double( 1.00000000000000000000e+00 ) +
			fX2 * (double( -1.66666671633720397949e-01 ) +
			fX2 * (double( 8.33333376795053482056e-03 ) +
			fX2 * (double( -1.98412497411482036114e-04 ) +
			fX2 * (double( 2.75565571428160183132e-06 ) +
			fX2 * (double( -2.50368472620721149724e-08 ) +
			fX2 * (double( 1.58849267073435385100e-10 ) +
			fX2 * double( -6.58925550841432672300e-13 ))))))));
}


constexpr float sin_adam42_s(float x)
{
    int32_t i32I = int32_t( x * (1.0f / PI) );
    x = (x - float( i32I ) * PI);

    float fX2 = x * x;
    float fX4 = fX2 * fX2;
    float fX6 = fX2 * fX4;
    float fX8 = fX4 * fX4;
    float fX10 = fX6 * fX4;
    float fX12 = fX6 * fX6;
    float fX14 = fX6 * fX8;

    return (i32I & 1) ?
        -x * (float( 1.00000000000000000000e+00 ) +
            (fX2 * float( -1.66666671633720397949e-01 )) +
            ((fX4 * float( 8.33333376795053482056e-03 )) +
            (fX6 * float( -1.98412497411482036114e-04 ))) +
            ((fX8 * float( 2.75565571428160183132e-06 )) +
            (fX10 * float( -2.50368472620721149724e-08 ))) +
            ((fX12 * float( 1.58849267073435385100e-10 )) +
            (fX14 * float( -6.58925550841432672300e-13 )))
                ):
        x * (float( 1.00000000000000000000e+00 ) +
            (fX2 * float( -1.66666671633720397949e-01 )) +
            ((fX4 * float( 8.33333376795053482056e-03 )) +
            (fX6 * float( -1.98412497411482036114e-04 ))) +
            ((fX8 * float( 2.75565571428160183132e-06 )) +
            (fX10 * float( -2.50368472620721149724e-08 ))) +
            ((fX12 * float( 1.58849267073435385100e-10 )) +
            (fX14 * float( -6.58925550841432672300e-13 )))
                );
}

constexpr double sin_adam42(double x)
{
    int32_t i32I = int32_t( x * (1.0f / PI) );
    x = (x - double( i32I ) * PI);

    double fX2 = x * x;
    double fX4 = fX2 * fX2;
    double fX6 = fX2 * fX4;
    double fX8 = fX4 * fX4;
    double fX10 = fX6 * fX4;
    double fX12 = fX6 * fX6;
    double fX14 = fX6 * fX8;

    return (i32I & 1) ?
        -x * (double( 1.00000000000000000000e+00 ) +
            (fX2 * double( -1.66666671633720397949e-01 )) +
            ((fX4 * double( 8.33333376795053482056e-03 )) +
            (fX6 * double( -1.98412497411482036114e-04 ))) +
            ((fX8 * double( 2.75565571428160183132e-06 )) +
            (fX10 * double( -2.50368472620721149724e-08 ))) +
            ((fX12 * double( 1.58849267073435385100e-10 )) +
            (fX14 * double( -6.58925550841432672300e-13 )))
                ):
        x * (double( 1.00000000000000000000e+00 ) +
            (fX2 * double( -1.66666671633720397949e-01 )) +
            ((fX4 * double( 8.33333376795053482056e-03 )) +
            (fX6 * double( -1.98412497411482036114e-04 ))) +
            ((fX8 * double( 2.75565571428160183132e-06 )) +
            (fX10 * double( -2.50368472620721149724e-08 ))) +
            ((fX12 * double( 1.58849267073435385100e-10 )) +
            (fX14 * double( -6.58925550841432672300e-13 )))
                );
}


double sin_spiro_c11( double _fX ) {
	int32_t i32I = int32_t( _fX * 0.31830988618379067153776752674503f );	// 1 / PI.
	_fX = (_fX - double( i32I ) * 3.1415926535897932384626433832795f);

	double fX2 = _fX * _fX;
	// Average error:
	//	0.000000071991304556528727213541666666667
	// Max error:
	//	0.000000543892383575439453125
	return (i32I & 1) ?
		-_fX * (double( 9.99999701976776123047e-01 ) +
			fX2 * (double( -1.66665777564048767090e-01 ) +
			fX2 * (double( 8.33255797624588012695e-03 ) +
			fX2 * (double( -1.98125766473822295666e-04 ) +
			fX2 * (double( 2.70405212177138309926e-06 ) +
			fX2 * double( -2.05329886426852681325e-08 )))))) :
		_fX * (double( 9.99999701976776123047e-01 ) +
			fX2 * (double( -1.66665777564048767090e-01 ) +
			fX2 * (double( 8.33255797624588012695e-03 ) +
			fX2 * (double( -1.98125766473822295666e-04 ) +
			fX2 * (double( 2.70405212177138309926e-06 ) +
			fX2 * double( -2.05329886426852681325e-08 ))))));
}
float sin_spiro_c11_s( float _fX ) {
	int32_t i32I = int32_t( _fX * 0.31830988618379067153776752674503f );	// 1 / PI.
	_fX = (_fX - float( i32I ) * 3.1415926535897932384626433832795f);

	float fX2 = _fX * _fX;
	// Average error:
	//	0.000000071991304556528727213541666666667
	// Max error:
	//	0.000000543892383575439453125
	return (i32I & 1) ?
		-_fX * (float( 9.99999701976776123047e-01 ) +
			fX2 * (float( -1.66665777564048767090e-01 ) +
			fX2 * (float( 8.33255797624588012695e-03 ) +
			fX2 * (float( -1.98125766473822295666e-04 ) +
			fX2 * (float( 2.70405212177138309926e-06 ) +
			fX2 * float( -2.05329886426852681325e-08 )))))) :
		_fX * (float( 9.99999701976776123047e-01 ) +
			fX2 * (float( -1.66665777564048767090e-01 ) +
			fX2 * (float( 8.33255797624588012695e-03 ) +
			fX2 * (float( -1.98125766473822295666e-04 ) +
			fX2 * (float( 2.70405212177138309926e-06 ) +
			fX2 * float( -2.05329886426852681325e-08 ))))));
}


class QPC
{
	uint64_t t1;
	uint64_t t2;
public:
	QPC() : t1(), t2() { }

	void start() { QueryPerformanceCounter((PLARGE_INTEGER) &t1); }
	void stop()  { QueryPerformanceCounter((PLARGE_INTEGER) &t2); }

	operator int() { return (int)(t2-t1); }
};


template<typename F> int test(F f, double low = -3.14, double high = +3.14)
{
	double sum = 0.0;
	double t = low;
	const double step = (high - low) / num_steps;

	QPC qpc;

	qpc.start();
		for(unsigned int i = num_steps; i --> 0; t += step)
			sum += f(t);
	qpc.stop();

	volatile double discard = sum;

	return qpc;
}


template<typename F> double func(F f, double t){ return std::abs(f(t)) - 1.0; }
template<typename T> int sign(T val) { return (T(0) < val) - (val < T(0)); }

template<typename F> double bisect(F f, double a, double b)
{
	constexpr unsigned int max_iter = 50;
	constexpr double tolerance = 0.00000001;

	unsigned int n = 0;
	double c;

	while(n++ < max_iter)
	{
		c = (a + b)/2;
		if(func(f, c) == 0.0 || std::abs(b-a) < tolerance)
			return c;
		if(sign(func(f,c)) == sign(func(f,a)))
			a = c;
		else
			b = c;
	}
	return c;
}


#define error(f) error_impl(#f, f)
template<typename F> void error_impl(const char* name, F f, double low = -3.14, double high = +3.14)
{
	printf("%-18s: ", name);
	double emax = 0.0;
	double eavg = 0.0;
	double sse  = 0.0;

	int    toobig  = 0;

	double t = low;

	const double step = (high - low) / num_steps;

	for(unsigned int i = num_steps; i --> 0; t += step)
	{
		double diff = std::abs(sin(t) - f(t));

		if(std::abs(f(t)) > 1.0) ++toobig;

		emax =  std::max(emax, diff);
		eavg += diff;
		sse  += diff*diff;
	}

	eavg /= num_steps;

	printf("emax=%1.15f eavg=%1.15f sse=%1.15f rmse=%1.15f   %d values > 1.0\n", emax, eavg, sse, sqrt(sse/num_steps), toobig);

	if(toobig)
	{
		double left, right;

		left  = bisect(f, -3*PI/2 - 0.1,    -3*PI/2);
		right = bisect(f, -3*PI/2,          -3*PI/2 + 0.1);
		printf("\tbracket around -3/2 pi where f(x) >  1 : [%+1.15f ; %+1.15f]    interval = %1.15f\n", left, right, std::abs(right - left));

		left  = bisect(f,   -PI/2 - 0.1,    -PI/2);
		right = bisect(f,   -PI/2,          -PI/2 + 0.1);
		printf("\tbracket around -pi/2   where f(x) < -1 : [%+1.15f ; %+1.15f]    interval = %1.15f\n", left, right, std::abs(right - left));

		left  = bisect(f,    PI/2 - 0.1,     PI/2);
		right = bisect(f,    PI/2,           PI/2 + 0.1);
		printf("\tbracket around  pi/2   where f(x) >  1 : [%+1.15f ; %+1.15f]    interval = %1.15f\n", left, right, std::abs(right - left));

		left  = bisect(f,  3*PI/2 - 0.1,     3*PI/2);
		right = bisect(f,  3*PI/2,           3*PI/2 + 0.1);
		printf("\tbracket around  3/2 pi where f(x) < -1 : [%+1.15f ; %+1.15f]    interval = %1.15f\n", left, right, std::abs(right - left));
	}
}


#define sanity(x) sanity_impl(#x, x)
template<typename F> void sanity_impl(const char* name, F f)
{
	{
		auto eq10 = [](double a, double b) { return std::abs(a -b) <= 0.00000000001; };

		auto err = [&name, &eq10](const char* msg, double expect, double value)
		{
			if(expect != value)
			{
				printf("%s %s failed the '%s' test   %s\n\texpect : %+1.16f\n\tvalue  : %+1.16f\n",
				       name, eq10(expect, value) ? "technically" : "", msg, eq10(expect, value) ? "(but has more than 10 good digits)": "", expect, value);
				printf("\tdiff   : %+1.16f\n",  std::abs(expect - value));
			}
		};

		if(f(-PI) != 0.0)    err("sin(-pi) == 0"   ,  0.0, f(-PI));
		if(f(-PI/2) != -1.0) err("sin(-pi/2) == -1", -1.0, f(-PI/2));
		if(f(0) != 0.0)      err("sin(0) == 0"     ,  0.0, f(0));
		if(f(PI/2) != 1.0)   err("sin(pi/2) == 1"  ,  1.0, f(PI/2));
		if(f(PI) != 0.0)     err("sin(pi) == 0"    ,  0.0, f(PI));
	}

	constexpr double eps = 0.001;
	uint64_t seed = 1234567890123ull;

	auto err = [&name, &eps, &f](const char* msg, double t) { printf("%s failed the '%s' test\n\tat t     = %1.16f\n\tf(t)     = %1.16f\n\tf(t+eps) = %1.16f\n", name, msg, t, f(t), f(t+eps)); };
	auto rnd = [](uint64_t x) { uint64_t z = (x += 0x9E3779B97F4A7C15ull); z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ull; z = (z ^ (z >> 27)) * 0x94D049BB133111EBull; return z ^ (z >> 31); };

	auto choose = [&eps, &seed, &rnd] (double low, double high) mutable { seed = rnd(seed); return low + ((seed & 1073741823)/1073741823.0 * (high - low - eps)); };

	int ctr = 0;
	constexpr int report_threshold = 3;

	auto incidents = [&ctr, &sani_steps]() mutable { if(ctr) printf("---> %'u incidents total (%'i samples)\n", ctr, sani_steps); ctr = 0; };

	for(int i = sani_steps; i --> 0;)
	{
		double t = choose(-2*PI, -3*PI/2);
		if(f(t + eps) <= f(t))
			if(++ctr <= report_threshold) err("monotonically rising in [-2 pi, -3/2 pi]", t);
	}
	incidents();
	for(int i = sani_steps; i --> 0;)
	{
		double t = choose(-3*PI/2, -PI/2);
		if(f(t + eps) >= f(t))
			if(++ctr <= report_threshold) err("monotonically falling in [-3/2 pi, -pi/2]", t);
	}
	incidents();
	for(int i = sani_steps; i --> 0;)
	{
		double t = choose(-PI/2, PI/2);
		if(f(t + eps) <= f(t))
			if(++ctr <= report_threshold) err("monotonically rising in [-pi/2, pi/2]", t);
	}
	incidents();
	for(int i = sani_steps; i --> 0;)
	{
		double t = choose(PI/2, 3*PI/2);
		if(f(t + eps) >= f(t))
			if(++ctr <= report_threshold) err("monotonically falling in [1/2 pi, 3/2 pi]", t);

	}
	incidents();
	for(int i = sani_steps; i --> 0;)
	{
		double t = choose(3*PI/2, 2*PI);
		if(f(t + eps) <= f(t))
			if(++ctr <= report_threshold) err("monotonically rising in [3/2 pi, 2 pi]", t);
	}
	incidents();
}


#include <cassert>
#include <locale.h>

int main()
{
	LARGE_INTEGER f;
	QueryPerformanceFrequency(&f);
	assert(f.HighPart == 0); // just to be sure

	setlocale(LC_ALL, "");
	printf("QPC resolution      = %'u ticks/sec\n", (unsigned int)f.LowPart);
	printf("Iterations per test = %'u\n", num_steps);

	{
		// avoid starting benchmark in low-power state
		volatile int wakeup_cpu = 0;
		for(unsigned int i = 0; i < 1'000'000; ++i)
			wakeup_cpu += (i*i) + (i-1)*i + (i+2)/(i+1);
	}

	int tsin           = test(sin);
	int tgarr11        = test(sin_garrett_c11);
	int tgarr11_s      = test(sin_garrett_c11_s);
	int tsin_spiro_s   = test(sin_spiro_s);
	int tsin_spiro     = test(sin_spiro);
	int tsin_adam42_s  = test(sin_adam42_s);
	int tsin_adam42    = test(sin_adam42);
	int tsin_spiro_c11_s = test(sin_spiro_c11_s);
	int tsin_spiro_c11   = test(sin_spiro_c11);

	puts("");
	auto show = [&qpc_per_sec = f.LowPart, &tsin](const char* name, int ticks)
	{
		printf("%-18s: %'-12d  [%f us per iteration] --> %-2.1f : 1\n", name, ticks, (((double) ticks) / num_steps  / qpc_per_sec) * 1'000'000.0, (double)tsin/ticks);
	};

	puts("\ntimings");
	show("sin",               tsin);
	show("sin_garrett_c11",   tgarr11);
	show("sin_garrett_c11_s", tgarr11_s);
	show("sin_spiro_s",       tsin_spiro_s);
	show("sin_spiro",         tsin_spiro);
	show("sin_adam42_s",      tsin_adam42_s);
	show("sin_adam42",        tsin_adam42);
	show("sin_spiro_c11_s",   tsin_spiro_c11_s);
	show("sin_spiro_c11",     tsin_spiro_c11);

	puts("\nerror metrics");

	error(sin);
	error(sin_garrett_c11);
	error(sin_garrett_c11_s);
	error(sin_spiro_s);
	error(sin_spiro);
	error(sin_adam42_s);
	error(sin_adam42);
	error(sin_spiro_c11_s);
	error(sin_spiro_c11);

	puts("\nsanity tests");
	sanity(sin);
	sanity(sin_garrett_c11);
	sanity(sin_garrett_c11_s);
	sanity(sin_spiro_s);
	sanity(sin_spiro);
	sanity(sin_adam42_s);
	sanity(sin_adam42);
	sanity(sin_spiro_c11_s);
	sanity(sin_spiro_c11);

    return 0;
}