Improved/formatted the hyperbolic and inverse trig functions

Author: ZERICO2005

src/libc/acoshf.c

@@ -1,8 +1,16 @@
 #include <math.h>
 
-float acoshf(float x)
-{
-    return logf(x + sqrtf(x * x - 1));
+/**
+ * @remarks Minimum ulp:
+ * ulp of -2  at +0x1.0007e8p+0 with ideal log1pf
+ * ulp of +10 at +0x1.01c964p+0 with current log1pf
+ */
+float acoshf(float x) {
+    if (x < 0x1.0p+63f) {
+        float t = x - 1.0f;
+        return log1pf(t + sqrtf(2.0f * t + t * t));
+    }
+    return logf(x) + (float)M_LN2;
 }
 
 double acosh(double) __attribute__((alias("acoshf")));

src/libc/acoshl.c

@@ -1,6 +1,28 @@
 #include <math.h>
+#include "__float64_constants.h"
 
-long double acoshl(long double x)
-{
-    return logl(x + sqrtl(x * x - 1));
+/* enable when log1pl is implemented */
+#if 0
+
+/**
+ * @remarks Minimum relative precision of:
+ * 2^-42.37 at +1.000000358e+00 with ideal log1pl
+ */
+long double acoshl(long double x) {
+    if (x < 0x1.0p+511L) {
+        long double t = x - 1.0L;
+        return log1pl(t + sqrtl(2.0L * t + t * t));
+    }
+    return logl(x) + F64_LN2;
 }
+
+#else
+
+long double acoshl(long double x) {
+    if (x < 0x1.0p+511L) {
+        return logl(x + sqrtl(x * x - 1.0L));
+    }
+    return logl(x) + F64_LN2;
+}
+
+#endif

src/libc/asinf.c

@@ -7,31 +7,33 @@
 
 #include <errno.h>
 #include <math.h>
+#include <stdbool.h>
 
-#define pio2  1.57079632679490
-
+/**
+ * @remarks Minimum ulp:
+ * ulp of +8 at +0x1.ffe956p-1
+ */
 float _asinf_c(float arg) {
-    float sign, temp;
+    bool arg_sign;
+    float temp;
+    arg_sign = signbit(arg);
+    arg = fabsf(arg);
 
-    sign = 1.;
-    if(arg < 0) {
-        arg = -arg;
-        sign = -1.;
-    }
-
-    if(arg > 1.) {
+    if(arg > 1.0f) {
         errno = EDOM;
-        return(0.);
+        return 0.0f;
     }
 
-    temp = sqrtf(1. - arg*arg);
-    if(arg > 0.7) {
-        temp = pio2 - atanf(temp/arg);
+    temp = sqrtf(1.0f - arg*arg);
+    if(arg > 0.7f) {
+        temp = (float)M_PI_2 - atanf(temp/arg);
     } else {
         temp = atanf(arg/temp);
     }
-
-    return(sign*temp);
+    if (arg_sign) {
+        temp = -temp;
+    }
+    return temp;
 }
 
 double _asin_c(double) __attribute__((alias("_asinf_c")));

src/libc/asinhf.c

@@ -1,8 +1,20 @@
 #include <math.h>
+#include "__float32_constants.h"
 
-float asinhf(float x)
-{
-    return logf(x + sqrtf(x * x + 1));
+/**
+ * @remarks Minimum ulp:
+ * ulp of -20 at +0x1.ff2afcp-5
+ */
+float asinhf(float arg) {
+    float x = fabsf(arg);
+    if (x < 0.0703125f) {
+        x = x - (x * x * x) * F32_1_DIV_6;
+    } else if (x < 0x1.0p+63f) {
+        x = logf(x + sqrtf(x * x + 1.0f));
+    } else {
+        x = logf(x) + F32_LN2;
+    }
+    return copysignf(x, arg);
 }
 
 double asinh(double) __attribute__((alias("asinhf")));

src/libc/asinhl.c

@@ -1,6 +1,18 @@
 #include <math.h>
+#include "__float64_constants.h"
 
-long double asinhl(long double x)
-{
-    return logl(x + sqrtl(x * x + 1));
+/**
+ * @remarks Relative precision of:
+ * 2^-42.41 at +9.768006857e-04 with ideal logl
+ */
+long double asinhl(long double arg) {
+    long double x = fabsl(arg);
+    if (x < 0x1.0p-10L) {
+        x = x - (x * x * x) * F64_1_DIV_6;
+    } else if (x < 0x1.0p+511L) {
+        x = logl(x + sqrtl(x * x + 1.0L));
+    } else {
+        x = logl(x) + F64_LN2;
+    }
+    return copysignl(x, arg);
 }

src/libc/atanhf.c

@@ -1,8 +1,15 @@
 #include <math.h>
 
-float atanhf(float x)
-{
-    return .5f * logf((1 + x) / (1 - x));
+/**
+ * @remarks Minimum ulp:
+ * ulp of +2  at +0x1.fe2dd0p-9 with ideal log1pf
+ * ulp of -9  at +0x1.e1e8d0p-5 with current log1pf and ideal logf
+ * ulp of +11 at +0x1.f30abcp-5 with current log1pf and logf
+ */
+float atanhf(float arg) {
+    float x = fabsf(arg);
+    x = 0.5f * log1pf((2.0f * x) / (1.0f - x));
+    return copysignf(x, arg);
 }
 
 double atanh(double) __attribute__((alias("atanhf")));

src/libc/atanhl.c

@@ -1,6 +1,23 @@
 #include <math.h>
 
+/* enable when log1pl is implemented */
+#if 0
+
+/**
+ * @remarks Relative precision of:
+ * 2^-51 at +1.277396381e-01 with ideal log1pl
+ */
+long double atanhl(long double arg) {
+    long double x = fabsl(arg);
+    x = 0.5L * log1pl((2.0L * x) / (1.0L - x));
+    return copysignl(x, arg);
+}
+
+#else
+
 long double atanhl(long double x)
 {
-    return .5L * logl((1 + x) / (1 - x));
+    return 0.5L * logl((1.0L + x) / (1.0L - x));
 }
+
+#endif

src/libc/coshf.c

@@ -22,30 +22,25 @@
 
 #include <math.h>
 
-#define p0  -0.630767364049772e+6
-#define p1  -0.899127202203951e+5
-#define p2  -0.289421135598956e+4
-#define p3  -0.263056321339750e+2
-#define q0  -0.630767364049772e+6
-#define q1  0.152151737879002e+5
-#define q2  -0.173678953558234e+3
-
+/**
+ * @remarks Minimum ulp:
+ * ulp of -1  at +0x1.39a3fep-12 with ideal expf (|x| < 21.0f)
+ * ulp of -18 at +0x1.0a049cp+4  with current expf (|x| < 21.0f)
+ */
 float _coshf_c(float arg) {
     float val;
-
-    if(arg < 0) {
-        arg = -arg;
-    }
+    arg = fabsf(arg);
 
     val = expf(arg);
 
-    if(arg > 21.) {
-        return val/2;
+    if (arg > 21.0f) {
+        return val / 2.0f;
     }
 
     val += expf(-arg);
-    val /= 2;
+    val /= 2.0f;
     return val;
 }
 
 double _cosh_c(double) __attribute__((alias("_coshf_c")));
+

src/libc/sinhf.c

@@ -22,40 +22,33 @@
 
 #include <math.h>
 
-#define p0  -0.630767364049772e+6
-#define p1  -0.899127202203951e+5
-#define p2  -0.289421135598956e+4
-#define p3  -0.263056321339750e+2
-#define q0  -0.630767364049772e+6
-#define q1  0.152151737879002e+5
-#define q2  -0.173678953558234e+3
-
+#define p0 -0.630767364049772e+6f
+#define p1 -0.899127202203951e+5f
+#define p2 -0.289421135598956e+4f
+#define p3 -0.263056321339750e+2f
+#define q0 -0.630767364049772e+6f
+#define q1  0.152151737879002e+5f
+#define q2 -0.173678953558234e+3f
+
+/**
+ * @remarks Minimum ulp:
+ * ulp of -3  at +0x1.eec25ap-9 with ideal expf (|x| < 21.0f)
+ * ulp of -18 at +0x1.0a049cp+4 with current expf (|x| < 21.0f)
+ */
 float _sinhf_c(float arg) {
-    float temp, argsq;
-    register int sign;
-
-    sign = 1;
-    if(arg < 0) {
-        arg = -arg;
-        sign = -1;
-    }
-
-    if(arg > 21.) {
-        temp = expf(arg)/2;
-        if (sign>0)
-            return(temp);
-        else
-            return(-temp);
+    float temp, argsq, x;
+    x = fabsf(arg);
+
+    if (x > 21.0f) {
+        temp = expf(x) / 2.0f;
+    } else if (x > 0.5f) {
+        temp = (expf(x) - expf(-x)) / 2.0f;
+    } else {
+        argsq = x * x;
+        temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0) * x;
+        temp /= (((argsq+q2)*argsq+q1)*argsq+q0);
     }
-
-    if(arg > 0.5) {
-        return(sign*(expf(arg) - expf(-arg))/2);
-    }
-
-    argsq = arg*arg;
-    temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg;
-    temp /= (((argsq+q2)*argsq+q1)*argsq+q0);
-    return(sign*temp);
+    return copysignf(temp, arg);
 }
 
 double _sinh_c(double, double *) __attribute__((alias("_sinhf_c")));