#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".
 * Corrected to avoid loss of accuracy for negative inputs.
 *
 * Max error: 1 ulp.
 * Avg error: 0.129 ulp.
 *
 * Works for all inputs.
 * Looks to be monotonic everywhere.
 *
 * [1] http://www.plunk.org/~hatch/rightway.php
 */
double plunk_atanh(double x)
{
	if (x == 1)
		return INFINITY;
	if (x == -1)
		return -INFINITY;
	if (x < 0)
		return -nr_log1p(-2*x/(1+x)) / 2;
	if (x > 0)
		return nr_log1p(2*x/(1-x)) / 2;
	return x;
}

/*
 * Algorithm using one round of Newton-Raphson.
 *
 * Max error: 2 ulp.
 * Avg error: 0.081 ulp.
 *
 * Works for all inputs.
 * Looks to be monotonic everywhere.
 */
double nr_atanh(double x)
{
	int sign;
	double y, t;

	if (x == 0)
		return x;
	if (x == 1)
		return INFINITY;
	if (x == -1)
		return -INFINITY;

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

	y = log((1+x)/(1-x)) / 2;
	t = tanh(y);

	return sign * (y - (t-x)/(1-t*t));
}

/*
 * Alternative algorithm using log1p in the obvious way.
 *
 * Max error: 1 ulp.
 * Avg error: 0.129 ulp.
 *
 * Works for all inputs.
 * Looks to be monotonic everywhere.
 * Nice and simple.
 */
double log1p_atanh(double x)
{
	if (x == 0)
		return x;
	if (x == 1)
		return INFINITY;
	if (x == -1)
		return -INFINITY;
	return (nr_log1p(x) - nr_log1p(-x)) / 2;
}

/* === TESTS === */

int num_tests = 0;
double max_plunk_err = 0.0;
double max_nr_err = 0.0;
double max_log1p_err = 0.0;
double total_plunk_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_atanh(double x)
{
	double y = atanh(x);
	double plunk_y = plunk_atanh(x);
	double nr_y = nr_atanh(x);
	double log1p_y = log1p_atanh(x);
	double plunk_err = err(y, plunk_y);
	double nr_err = err(y, nr_y);
	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 > 1 || (fn)(next_x) >= (fn)(x) ) ? \
	"Monotonic" : "Not monotonic"

	printf("x = %.20g\n", x);
	printf("Standard result: %.20g %s\n", y, monotonic(atanh));
	printf("Plunk atanh():   %.20g,  err=%f %s\n", plunk_y, plunk_err,
		   monotonic(plunk_atanh));
	printf("NR atanh():      %.20g,  err=%f %s\n", nr_y, nr_err,
		   monotonic(nr_atanh));
	printf("log1p atanh():   %.20g,  err=%f %s\n", log1p_y, log1p_err,
		   monotonic(log1p_atanh));

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

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

int main()
{
	int s;

	for (s=-1; s<=1; s+=2)
	{
		test_atanh(s*0.0);
		test_atanh(s*nextafter(0, 1));
		test_atanh(s*4.94e-324);
		test_atanh(s*9.9e-324);
		test_atanh(s*2.23e-320);
		test_atanh(s*2.23e-315);
		test_atanh(s*2.23e-308);
		test_atanh(s*3.141e-300);
		test_atanh(s*3.141e-200);
		test_atanh(s*3.141e-100);
		test_atanh(s*3.141e-90);
		test_atanh(s*3.141e-80);
		test_atanh(s*3.141e-70);
		test_atanh(s*3.141e-60);
		test_atanh(s*3.141e-50);
		test_atanh(s*3.141e-40);
		test_atanh(s*3.141e-30);
		test_atanh(s*3.141e-20);
		test_atanh(s*3.141e-19);
		test_atanh(s*3.141e-18);
		test_atanh(s*3.141e-17);
		test_atanh(s*3.141e-16);
		test_atanh(s*3.141e-15);
		test_atanh(s*3.141e-14);
		test_atanh(s*3.141e-13);
		test_atanh(s*3.141e-12);
		test_atanh(s*3.141e-11);
		test_atanh(s*3.141e-10);
		test_atanh(s*3.141e-9);
		test_atanh(s*3.141e-8);
		test_atanh(s*3.141e-7);
		test_atanh(s*3.141e-6);
		test_atanh(s*3.141e-5);
		test_atanh(s*3.141e-4);
		test_atanh(s*3.141e-3);
		test_atanh(s*3.141e-2);
		test_atanh(s*0.13141);
		test_atanh(s*0.23141);
		test_atanh(s*0.33141);
		test_atanh(s*0.43141);
		test_atanh(s*0.53141);
		test_atanh(s*0.63141);
		test_atanh(s*0.73141);
		test_atanh(s*0.83141);
		test_atanh(s*0.93141);
		test_atanh(s*(1-3.141e-2));
		test_atanh(s*(1-3.141e-3));
		test_atanh(s*(1-3.141e-4));
		test_atanh(s*(1-3.141e-5));
		test_atanh(s*(1-3.141e-6));
		test_atanh(s*(1-3.141e-7));
		test_atanh(s*(1-3.141e-8));
		test_atanh(s*(1-3.141e-9));
		test_atanh(s*(1-3.141e-10));
		test_atanh(s*(1-3.141e-11));
		test_atanh(s*(1-3.141e-12));
		test_atanh(s*(1-3.141e-13));
		test_atanh(s*(1-3.141e-14));
		test_atanh(s*(1-3.141e-15));
		test_atanh(s*(1-3.141e-16));
		test_atanh(s*nextafter(1, 0));
		test_atanh(s);
	}

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

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

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

	return 0;
}