#include <math.h>
#include <stdio.h>
#include <string.h>
#include <sys/time.h>

#define REPS 100000000

/* mercators */

inline double merc1(double x) {
	return asinh(tan(x));
}

inline double merc2(double x) {
	return atanh(sin(x));
}

inline double merc3(double x) {
	return 2.0*atanh(tan(x/2.0));
}

inline double merc4(double x) {
	return log(tan(M_PI/4.0 + x/2.0));
}

inline double merc5(double x) {
	return 0.5*log((1.0+sin(x))/(1.0-sin(x)));
}

inline double merc6(double x) {
	return log(tan(x) + 1.0/cos(x));
}

inline double merc7(double x) {
	return log((1.0+sin(x))/cos(x));
}

/* reverse mercators */
inline double revmerc1(double x) {
	return atan(sinh(x));
}

inline double revmerc2(double x) {
	return asin(tanh(x));
}

inline double revmerc3(double x) {
	return 2.0*atan(tanh(0.5*x));
}

inline double revmerc4(double x) {
	return 2.0*atan(exp(x)) - M_PI/2.0;
}

#define FUNCTEST(function, reverse, testvals) \
	{ \
		int i; \
		struct timeval start, end; \
		double sum = 0.0, error = 0.0; \
		gettimeofday(&start, 0); \
		for (i = 0; i < REPS+1; i++) \
			sum += function((double)i/(double)REPS-0.5); \
		gettimeofday(&end, 0); \
		for (i = 0; i < sizeof(testvals)/sizeof(testvals[0]); i++) \
			error += reverse(function(testvals[i])) - testvals[i]; \
		printf("%s: %f seconds, f(0.5)=%f, sum=%le, error=%le\n", #function, (float)(end.tv_sec - start.tv_sec) + (float)(end.tv_usec - start.tv_usec)/1000000.0f, function(0.5), fabs(sum), fabs(error)); \
	}

int matrix() {
	double (*mercs[])(double) = { merc1, merc2, merc3, merc4, merc5, merc6, merc7 };
	double (*revmercs[])(double) = { revmerc1, revmerc2, revmerc3, revmerc4 };

	double mercerrs[sizeof(mercs)/sizeof(mercs[0])] = {0,};
	double revmercerrs[sizeof(mercs)/sizeof(mercs[0])] = {0,};

	int i,j;
	double a;

	printf("     ");
	for (j = 0; j < sizeof(revmercs)/sizeof(revmercs[0]); j++)
		printf("      revmerc%d", j+1);
	printf("\n");

	for (i = 0; i < sizeof(mercs)/sizeof(mercs[0]); i++) {
		printf("merc%d", i+1);
		for (j = 0; j < sizeof(revmercs)/sizeof(revmercs[0]); j++) {
			double error = 0.0;
			for (a = -1.57; a < 1.57; a += 0.001)
				error += fabs(revmercs[j](mercs[i](a)) - a);
			printf("%14e", error);
			mercerrs[i] += error;
			revmercerrs[j] += error;
		}
		printf("\n");
	}

	printf("\n");

	for (i = 0; i < sizeof(mercs)/sizeof(mercs[0]); i++)
		printf("merc%d: %e\n", i+1, mercerrs[i]);

	printf("\n");

	for (i = 0; i < sizeof(revmercs)/sizeof(revmercs[0]); i++)
		printf("revmerc%d: %e\n", i+1, revmercerrs[i]);

	return 0;
}

int benchmark() {
	double angle_vals[] = { -1.57, -1.0, 0.0, 1.0, 1.57 };
	double merc_vals[] = { -1.0, -0.1, 0.0, 0.1, 1.0 };

	/* mercator */
	FUNCTEST(merc1, revmerc1, angle_vals);
	FUNCTEST(merc2, revmerc2, angle_vals);
	FUNCTEST(merc3, revmerc3, angle_vals);
	FUNCTEST(merc4, revmerc4, angle_vals);
	FUNCTEST(merc5, revmerc4, angle_vals);
	FUNCTEST(merc6, revmerc4, angle_vals);
	FUNCTEST(merc7, revmerc4, angle_vals);

	/* reverse mercator */
	FUNCTEST(revmerc1, merc1, merc_vals);
	FUNCTEST(revmerc2, merc2, merc_vals);
	FUNCTEST(revmerc3, merc3, merc_vals);
	FUNCTEST(revmerc4, merc4, merc_vals);

	return 0;
}

int main(int argc, char **argv) {
	if (argc == 2 && strcmp(argv[1], "-m") == 0)
		return matrix();
	else if (argc == 2 && strcmp(argv[1], "-b") == 0)
		return benchmark();

	fprintf(stderr, "Usage: %s -m|-b\n", argv[0]);
	return 1;
}