On Sun, 3 Feb 2019 at 15:12, Tom Lane <[email protected]> wrote:
>
> Andrew Gierth <[email protected]> writes:
> > The spec doesn't require the inverse functions (asinh, acosh, atanh),
> > but surely there is no principled reason to omit them?
>
> +1 --- AFAICS, the C library has offered all six since C89.
>
+1 for including the inverse functions. However, it looks to me like
the inverse functions are C99-specific, so they might not be available
on all supported platforms. If they're not, we may need to provide our
own implementations.
I did a bit of research and had play. It looks like it might not be
too hard to provide our own implementations, but it does take a bit of
care to avoid rounding and overflow errors. Attached are some
standalone C programs where I tested various algorithms. A decent
approach seems to be to either use log1p() (which is itself
C99-specific, hence that's the first thing I played with) or to use a
single round of Newton-Raphson to correct rounding errors, in a
similar way to how we implement cbrt() on platforms that don't have
that.
Of course that may all be moot -- those functions may in fact be
available everywhere we care about, but it was interesting to play
around with them anyway.
Regards,
Dean
#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;
}
#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;
}
#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;
}
#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;
}
