#include <float.h>
#include <math.h>
#include <stdio.h>

double nr_log1p(double x)
{
	double y, e;

	if (x == 0)
		return x;

	y = log(1+x);
	e = exp(y);

	return y - (e-1-x)/e;
}

/*
 * Algorithm recommended by [1] "The Right Way to Calculate Stuff".
 *
 * Loses all accuracy for negative inputs, but easily fixed, as below.
 * Also loses all accuracy for large inputs. Not so easily fixed.
 *
 * Max error: 12193974156573 ulp (when x ~ +/- 5.6e200).
 * Avg error: 554271552571 ulp.
 *
 * Not the right way to calculate asinh.
 *
 * [1] http://www.plunk.org/~hatch/rightway.php
 */
double plunk_asinh(double x)
{
	if (x < 0)
		return -nr_log1p(-x * (1 - x/(sqrt(x*x+1)+1)));
	if (x > 0)
		return nr_log1p(x * (1 + x/(sqrt(x*x+1)+1)));
	return x;
}

/*
 * Alternative algorithm using one round of Newton-Raphson.
 * Have to stitch together 2 different formulae at x=1 to avoid overflow risk.
 *
 * Max error: 1 ulp.
 * Avg error: 0.117 ulp.
 *
 * Works for all inputs.
 * Looks to be monotonic everywhere.
 */
double nr_asinh(double x)
{
	int sign;
	double y;

	if (x == 0)
		return x;

	sign = x < 0 ? -1 : 1;
	x = fabs(x);

	if (x < 1)
		y = log(x+sqrt(1+x*x));
	else
		y = log(x) + log(1+sqrt(1+1/(x*x)));

	return sign * (y - (sinh(y)-x)/cosh(y));
}

/* === TESTS === */

int num_tests = 0;
double max_plunk_err = 0.0;
double max_nr_err = 0.0;
double total_plunk_err = 0.0;
double total_nr_err = 0.0;

double err(double x, double y)
{
	if (x < y)
	{
		double diff = y - x;
		double ulp = nextafter(x, y) - x;
		return diff / ulp;
	}
	if (x > y)
	{
		double diff = x - y;
		double ulp = nextafter(y, x) - y;
		return diff / ulp;
	}
	return 0.0;
}

void test_asinh(double x)
{
	double y = asinh(x);
	double plunk_y = plunk_asinh(x);
	double nr_y = nr_asinh(x);
	double plunk_err = err(y, plunk_y);
	double nr_err = err(y, nr_y);

	double prev_x = nextafter(x, -DBL_MAX);
	double next_x = nextafter(x, DBL_MAX);
#define monotonic(fn) \
	( prev_x == -DBL_MAX || (fn)(prev_x) <= (fn)(x) ) && \
	( next_x == DBL_MAX || (fn)(next_x) >= (fn)(x) ) ? \
	"Monotonic" : "Not monotonic"

	printf("x = %.20g\n", x);
	printf("Standard result: %.20g %s\n", y, monotonic(asinh));
	printf("Plunk asinh():   %.20g,  err=%f %s\n", plunk_y, plunk_err,
		   monotonic(plunk_asinh));
	printf("NR asinh():      %.20g,  err=%f %s\n", nr_y, nr_err,
		   monotonic(nr_asinh));

	if (plunk_err > max_plunk_err) max_plunk_err = plunk_err;
	if (nr_err > max_nr_err) max_nr_err = nr_err;

	total_plunk_err += plunk_err;
	total_nr_err += nr_err;
	num_tests++;
}

int main()
{
	int s;

	for (s=-1; s<=1; s+=2)
	{
		test_asinh(s*0.0);
		test_asinh(s*nextafter(0, 1));
		test_asinh(s*4.94e-324);
		test_asinh(s*9.9e-324);
		test_asinh(s*2.23e-320);
		test_asinh(s*2.23e-315);
		test_asinh(s*2.23e-308);
		test_asinh(s*3.141e-300);
		test_asinh(s*3.141e-200);
		test_asinh(s*3.141e-100);
		test_asinh(s*3.141e-90);
		test_asinh(s*3.141e-80);
		test_asinh(s*3.141e-70);
		test_asinh(s*3.141e-60);
		test_asinh(s*3.141e-50);
		test_asinh(s*3.141e-40);
		test_asinh(s*3.141e-30);
		test_asinh(s*3.141e-20);
		test_asinh(s*3.141e-19);
		test_asinh(s*3.141e-18);
		test_asinh(s*3.141e-17);
		test_asinh(s*3.141e-16);
		test_asinh(s*3.141e-15);
		test_asinh(s*3.141e-14);
		test_asinh(s*3.141e-13);
		test_asinh(s*3.141e-12);
		test_asinh(s*3.141e-11);
		test_asinh(s*3.141e-10);
		test_asinh(s*3.141e-9);
		test_asinh(s*3.141e-8);
		test_asinh(s*3.141e-7);
		test_asinh(s*3.141e-6);
		test_asinh(s*3.141e-5);
		test_asinh(s*3.141e-4);
		test_asinh(s*3.141e-3);
		test_asinh(s*3.141e-2);
		test_asinh(s*0.3141);
		test_asinh(s*0.6324432);
		test_asinh(s*0.96324432);
		test_asinh(s*nextafter(1, 0));
		test_asinh(s*1.0);
		test_asinh(s*nextafter(1, 2));
		test_asinh(s*1.6324432);
		test_asinh(s*2.6324432);
		test_asinh(s*5.6324432);
		test_asinh(s*24.6324432);
		test_asinh(s*5.6324432e1);
		test_asinh(s*5.6324432e2);
		test_asinh(s*5.6324432e3);
		test_asinh(s*5.6324432e4);
		test_asinh(s*5.6324432e5);
		test_asinh(s*5.6324432e6);
		test_asinh(s*5.6324432e7);
		test_asinh(s*5.6324432e8);
		test_asinh(s*5.6324432e9);
		test_asinh(s*5.6324432e10);
		test_asinh(s*5.6324432e11);
		test_asinh(s*5.6324432e12);
		test_asinh(s*5.6324432e13);
		test_asinh(s*5.6324432e14);
		test_asinh(s*5.6324432e15);
		test_asinh(s*5.6324432e16);
		test_asinh(s*5.6324432e17);
		test_asinh(s*5.6324432e18);
		test_asinh(s*5.6324432e19);
		test_asinh(s*5.6324432e20);
		test_asinh(s*5.6324432e30);
		test_asinh(s*5.6324432e40);
		test_asinh(s*5.6324432e50);
		test_asinh(s*5.6324432e60);
		test_asinh(s*5.6324432e70);
		test_asinh(s*5.6324432e80);
		test_asinh(s*5.6324432e90);
		test_asinh(s*5.6324432e100);
		test_asinh(s*5.6324432e200);
		test_asinh(s*5.6324432e300);
		test_asinh(s*5.6324432e307);
		test_asinh(s*8.98e307);
		test_asinh(s*1.797e308);
		test_asinh(s*DBL_MAX);
	}

	printf("\n");
	printf("Number of tests: %d\n", num_tests);

	printf("\n");
	printf("Max error for Plunk asinh(): %f\n", max_plunk_err);
	printf("Max error for NR asinh():    %f\n", max_nr_err);

	printf("\n");
	printf("Avg. error for Plunk asinh(): %f\n", total_plunk_err / num_tests);
	printf("Avg. error for NR asinh():    %f\n", total_nr_err / num_tests);

	return 0;
}