#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 using one round of Newton-Raphson.
 * Not expected to work well, because of infinite slope at x=1.
 *
 * Max error: 25216051 ulp (near x=1).
 * Avg error: 1236022 ulp.
 *
 * Not recommended.
 */
double nr_acosh(double x)
{
	double y;

	if (x < 1)
		return NAN;
	if (x == 1)
        return 0;

	y = log(x) + log(1+sqrt(1-1/(x*x)));
	return y - (cosh(y)-x)/sinh(y);
}

/*
 * Alternative algorithm using log1p.
 * Need to be careful to avoid overflow risk for large inputs.
 *
 * Max error: 2 ulp.
 * Avg error: 0.210 ulp.
 *
 * Works for all inputs.
 * Looks to be monotonic everywhere.
 */
double log1p_acosh(double x)
{
	double t, u;

	if (x < 1)
		return NAN;
	if (x == 1)
		return 0;

	t = x-1;
	u = t/x;
	return log(x) + nr_log1p(u * sqrt(1+2/t));
}

/* === TESTS === */

int num_tests = 0;
double max_nr_err = 0.0;
double max_log1p_err = 0.0;
double total_nr_err = 0.0;
double total_log1p_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_acosh(volatile double x)
{
	volatile double y = acosh(x);
	volatile double nr_y = nr_acosh(x);
	volatile double log1p_y = log1p_acosh(x);
	volatile double nr_err = err(y, nr_y);
	volatile double log1p_err = err(y, log1p_y);

	double prev_x = nextafter(x, -DBL_MAX);
	double next_x = nextafter(x, DBL_MAX);
#define monotonic(fn) \
	( prev_x < 1 || (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(acosh));
	printf("NR acosh():      %.20g,  err=%f %s\n", nr_y, nr_err,
		   monotonic(nr_acosh));
	printf("log1p acosh():   %.20g,  err=%f %s\n", log1p_y, log1p_err,
		   monotonic(log1p_acosh));

	if (nr_err > max_nr_err) max_nr_err = nr_err;
	if (log1p_err > max_log1p_err) max_log1p_err = log1p_err;

	total_nr_err += nr_err;
	total_log1p_err += log1p_err;
	num_tests++;
}

int main()
{
	test_acosh(1);
	test_acosh(1+3.141e-20);
	test_acosh(1+3.141e-19);
	test_acosh(1+3.141e-18);
	test_acosh(1+3.141e-17);
	test_acosh(1+3.141e-16);
	test_acosh(nextafter(1, 2));
	test_acosh(1+3.141e-15);
	test_acosh(1+3.141e-14);
	test_acosh(1+3.141e-13);
	test_acosh(1+3.141e-12);
	test_acosh(1+3.141e-11);
	test_acosh(1+3.141e-10);
	test_acosh(1+3.141e-9);
	test_acosh(1+3.141e-8);
	test_acosh(1+3.141e-7);
	test_acosh(1+3.141e-6);
	test_acosh(1+3.141e-5);
	test_acosh(1+3.141e-4);
	test_acosh(1+3.141e-3);
	test_acosh(1+3.141e-2);
	test_acosh(1.3141);
	test_acosh(1.6324432);
	test_acosh(1.96324432);
	test_acosh(2.6324432);
	test_acosh(3.6324432);
	test_acosh(5.6324432);
	test_acosh(24.6324432);
	test_acosh(5.6324432e1);
	test_acosh(5.6324432e2);
	test_acosh(5.6324432e3);
	test_acosh(5.6324432e4);
	test_acosh(5.6324432e5);
	test_acosh(5.6324432e6);
	test_acosh(5.6324432e7);
	test_acosh(5.6324432e8);
	test_acosh(5.6324432e9);
	test_acosh(5.6324432e10);
	test_acosh(5.6324432e11);
	test_acosh(5.6324432e12);
	test_acosh(5.6324432e13);
	test_acosh(5.6324432e14);
	test_acosh(5.6324432e15);
	test_acosh(5.6324432e16);
	test_acosh(5.6324432e17);
	test_acosh(5.6324432e18);
	test_acosh(5.6324432e19);
	test_acosh(5.6324432e20);
	test_acosh(5.6324432e30);
	test_acosh(5.6324432e40);
	test_acosh(5.6324432e50);
	test_acosh(5.6324432e60);
	test_acosh(5.6324432e70);
	test_acosh(5.6324432e80);
	test_acosh(5.6324432e90);
	test_acosh(5.6324432e100);
	test_acosh(5.6324432e200);
	test_acosh(5.6324432e300);
	test_acosh(5.6324432e307);
	test_acosh(8.98e307);
	test_acosh(1.797e308);
	test_acosh(DBL_MAX);

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

	printf("\n");
	printf("Max error for NR acosh():    %f\n", max_nr_err);
	printf("Max error for log1p acosh(): %f\n", max_log1p_err);

	printf("\n");
	printf("Avg. error for NR acosh():    %f\n", total_nr_err / num_tests);
	printf("Avg. error for log1p acosh(): %f\n", total_log1p_err / num_tests);

	return 0;
}