math: rewrite inverse hyperbolic functions to be simpler/smaller

[musl] / src / math / asinh.c

diff --git a/src/math/asinh.c b/src/math/asinh.c
index 11bbd71..4152dc3 100644 (file)
--- a/src/math/asinh.c
+++ b/src/math/asinh.c
@@ -1,55 +1,28 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_asinh.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.
- * ====================================================
- */
-/* asinh(x)
- * Method :
- *      Based on
- *              asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
- *      we have
- *      asinh(x) := x  if  1+x*x=1,
- *               := sign(x)*(log(x)+ln2)) for large |x|, else
- *               := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
- *               := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
- */
-
 #include "libm.h"
 
-static const double
-ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
-huge= 1.00000000000000000000e+300;
-
+/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
 double asinh(double x)
 {
-       double t,w;
-       int32_t hx,ix;
+       union {double f; uint64_t i;} u = {.f = x};
+       unsigned e = u.i >> 52 & 0x7ff;
+       unsigned s = u.i >> 63;
 
-       GET_HIGH_WORD(hx, x);
-       ix = hx & 0x7fffffff;
-       if (ix >= 0x7ff00000)   /* x is inf or NaN */
-               return x+x;
-       if (ix < 0x3e300000) {  /* |x| < 2**-28 */
-               /* return x inexact except 0 */
-               if (huge+x > 1.0)
-                       return x;
-       }
-       if (ix > 0x41b00000) {  /* |x| > 2**28 */
-               w = log(fabs(x)) + ln2;
-       } else if (ix > 0x40000000) {  /* 2**28 > |x| > 2.0 */
-               t = fabs(x);
-               w = log(2.0*t + 1.0/(sqrt(x*x+1.0)+t));
-       } else {                /* 2.0 > |x| > 2**-28 */
-               t = x*x;
-               w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t)));
+       /* |x| */
+       u.i &= (uint64_t)-1/2;
+       x = u.f;
+
+       if (e >= 0x3ff + 26) {
+               /* |x| >= 0x1p26 or inf or nan */
+               x = log(x) + 0.693147180559945309417232121458176568;
+       } else if (e >= 0x3ff + 1) {
+               /* |x| >= 2 */
+               x = log(2*x + 1/(sqrt(x*x+1)+x));
+       } else if (e >= 0x3ff - 26) {
+               /* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */
+               x = log1p(x + x*x/(sqrt(x*x+1)+1));
+       } else {
+               /* |x| < 0x1p-26, raise inexact if x != 0 */
+               FORCE_EVAL(x + 0x1p1000);
        }
-       if (hx > 0)
-               return w;
-       return -w;
+       return s ? -x : x;
 }