math: rewrite inverse hyperbolic functions to be simpler/smaller

modifications:
* avoid unsigned->signed integer conversion
* do not handle special cases when they work correctly anyway
* more strict threshold values (0x1p26 instead of 0x1p28 etc)
* smaller code, cleaner branching logic
* same precision as the old code:
    acosh(x) has up to 2ulp error in [1,1.125]
    asinh(x) has up to 1.6ulp error in [0.125,0.5], [-0.5,-0.125]
    atanh(x) has up to 1.7ulp error in [0.125,0.5], [-0.5,-0.125]

1 files changed, 12 insertions, 46 deletions

diff --git a/src/math/acoshl.c b/src/math/acoshl.c
index a4024516..472c71cb 100644
--- a/src/math/acoshl.c
+++ b/src/math/acoshl.c
@@ -1,28 +1,3 @@
-/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_acoshl.c */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* acoshl(x)
- * Method :
- *      Based on
- *              acoshl(x) = logl [ x + sqrtl(x*x-1) ]
- *      we have
- *              acoshl(x) := logl(x)+ln2,       if x is large; else
- *              acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
- *              acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
- *
- * Special cases:
- *      acoshl(x) is NaN with signal if x<1.
- *      acoshl(NaN) is NaN without signal.
- */
-
 #include "libm.h"
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
@@ -31,29 +6,20 @@ long double acoshl(long double x)
 	return acosh(x);
 }
 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double
-ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
-
+/* acosh(x) = log(x + sqrt(x*x-1)) */
 long double acoshl(long double x)
 {
-	long double t;
-	uint32_t se,i0,i1;
+	union {
+		long double f;
+		struct{uint64_t m; int16_t se; uint16_t pad;} i;
+	} u = {.f = x};
 
-	GET_LDOUBLE_WORDS(se, i0, i1, x);
-	if (se < 0x3fff || se & 0x8000) {  /* x < 1 */
-		return (x-x)/(x-x);
-	} else if (se >= 0x401d) {  /* x > 2**30 */
-		if (se >= 0x7fff)  /* x is inf or NaN */
-			return x+x;
-		return logl(x) + ln2;  /* acoshl(huge) = logl(2x) */
-	} else if (((se-0x3fff)|i0|i1) == 0) {
-		return 0.0;            /* acosh(1) = 0 */
-	} else if (se > 0x4000) {  /* x > 2 */
-		t = x*x;
-		return logl(2.0*x - 1.0/(x + sqrtl(t - 1.0)));
-	}
-	/* 1 < x <= 2 */
-	t = x - 1.0;
-	return log1pl(t + sqrtl(2.0*t + t*t));
+	if (u.i.se < 0x3fff + 1)
+		/* x < 2, invalid if x < 1 or nan */
+		return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1)));
+	if (u.i.se < 0x3fff + 32)
+		/* x < 0x1p32 */
+		return logl(2*x - 1/(x+sqrtl(x*x-1)));
+	return logl(x) + 0.693147180559945309417232121458176568L;
 }
 #endif