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

/*
 * A commonly used algorithm, due to Kahan (citation needed).
 *
 * Max error: 1 ulp.
 * Avg error: 0.158 ulp.
 *
 * It's not obvious, but this does appear to be monotonic at the cutover point
 * (at least on my machine). Can that be relied upon on all platforms?
 *
 * -5.5511151231257851673e-17 -> -5.5511151231257851673e-17 (u != 1 branch)
 * -5.5511151231257839347e-17 -> -5.5511151231257839347e-17 (u != 1 branch)
 * -5.5511151231257827021e-17 -> -5.5511151231257827021e-17 (u == 1 branch)
 * -5.5511151231257820858e-17 -> -5.5511151231257820858e-17 (u == 1 branch)
 *
 * 1.1102230246251564172e-16 -> 1.1102230246251564172e-16 (u == 1 branch)
 * 1.1102230246251565404e-16 -> 1.1102230246251565404e-16 (u == 1 branch)
 * 1.1102230246251567869e-16 -> 1.1102230246251565404e-16 (u != 1 branch)
 * 1.1102230246251570335e-16 -> 1.1102230246251567869e-16 (u != 1 branch)
 *
 * However, it doesn't appear to be monotonic at various other points.
 */
double kahan_log1p(double x)
{
	volatile double u = 1+x;
	return u == 1 ? x : x * (log(u) / (u-1));
}

/*
 * Algorithm used in the GNU Scientific Library.
 *
 * Max error: 1 ulp.
 * Avg error: 0.086 ulp.
 *
 * Where does this formula come from? Licensing?
 * Doesn't return -0 when x is -0, but that could be fixed up.
 * Looks to be monotonic everywhere.
 */
double gsl_log1p(double x)
{
	volatile double y, z;
	y = 1 + x;
	z = y - 1;
	return log(y) - (z-x)/y;
}

/*
 * Alternative algorithm using one round of Newton-Raphson.
 *
 * Max error: 1 ulp.
 * Avg error: 0.094 ulp.
 *
 * Works well for all inputs.
 * Looks to be monotonic everywhere.
 */
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;
}

/* === TESTS === */

int num_tests = 0;
double max_kahan_err = 0.0;
double max_gsl_err = 0.0;
double max_nr_err = 0.0;
double total_kahan_err = 0.0;
double total_gsl_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_log1p(double x)
{
	double y = log1p(x);
	double kahan_y = kahan_log1p(x);
	double gsl_y = gsl_log1p(x);
	double nr_y = nr_log1p(x);
	double kahan_err = err(y, kahan_y);
	double gsl_err = err(y, gsl_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 == -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(log1p));
	printf("Kahan log1p():   %.20g,  err=%f %s\n", kahan_y, kahan_err,
		   monotonic(kahan_log1p));
	printf("GSL log1p():     %.20g,  err=%f %s\n", gsl_y, gsl_err,
		   monotonic(gsl_log1p));
	printf("NR log1p():      %.20g,  err=%f %s\n", nr_y, nr_err,
		   monotonic(nr_log1p));

	if (kahan_err > max_kahan_err) max_kahan_err = kahan_err;
	if (gsl_err > max_gsl_err) max_gsl_err = gsl_err;
	if (nr_err > max_nr_err) max_nr_err = nr_err;

	total_kahan_err += kahan_err;
	total_gsl_err += gsl_err;
	total_nr_err += nr_err;
	num_tests++;
}

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

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

	printf("\n");
	printf("Max error for Kahan log1p(): %f\n", max_kahan_err);
	printf("Max error for GSL log1p():   %f\n", max_gsl_err);
	printf("Max error for NR log1p():    %f\n", max_nr_err);

	printf("\n");
	printf("Avg. error for Kahan log1p(): %f\n", total_kahan_err / num_tests);
	printf("Avg. error for GSL log1p():   %f\n", total_gsl_err / num_tests);
	printf("Avg. error for NR log1p():    %f\n", total_nr_err / num_tests);

	return 0;
}