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, 21 insertions, 41 deletions

diff --git a/src/math/asinhf.c b/src/math/asinhf.c
index efe3af94..fc9f0911 100644
--- a/src/math/asinhf.c
+++ b/src/math/asinhf.c
@@ -1,48 +1,28 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_asinhf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
 #include "libm.h"
 
-static const float
-ln2 = 6.9314718246e-01, /* 0x3f317218 */
-huge= 1.0000000000e+30;
-
+/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
 float asinhf(float x)
 {
-	float t,w;
-	int32_t hx,ix;
+	union {float f; uint32_t i;} u = {.f = x};
+	uint32_t i = u.i & 0x7fffffff;
+	unsigned s = u.i >> 31;
 
-	GET_FLOAT_WORD(hx, x);
-	ix = hx & 0x7fffffff;
-	if (ix >= 0x7f800000)   /* x is inf or NaN */
-		return x+x;
-	if (ix < 0x31800000) {  /* |x| < 2**-28 */
-		/* return x inexact except 0 */
-		if (huge+x > 1.0f)
-			return x;
-	}
-	if (ix > 0x4d800000) {  /* |x| > 2**28 */
-		w = logf(fabsf(x)) + ln2;
-	} else if (ix > 0x40000000) {  /* 2**28 > |x| > 2.0 */
-		t = fabsf(x);
-		w = logf(2.0f*t + 1.0f/(sqrtf(x*x+1.0f)+t));
-	} else {                /* 2.0 > |x| > 2**-28 */
-		t = x*x;
-		w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t)));
+	/* |x| */
+	u.i = i;
+	x = u.f;
+
+	if (i >= 0x3f800000 + (12<<23)) {
+		/* |x| >= 0x1p12 or inf or nan */
+		x = logf(x) + 0.693147180559945309417232121458176568f;
+	} else if (i >= 0x3f800000 + (1<<23)) {
+		/* |x| >= 2 */
+		x = logf(2*x + 1/(sqrtf(x*x+1)+x));
+	} else if (i >= 0x3f800000 - (12<<23)) {
+		/* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */
+		x = log1pf(x + x*x/(sqrtf(x*x+1)+1));
+	} else {
+		/* |x| < 0x1p-12, raise inexact if x!=0 */
+		FORCE_EVAL(x + 0x1p120f);
 	}
-	if (hx > 0)
-		return w;
-	return -w;
+	return s ? -x : x;
 }