Index: fusl/src/math/sinh.c |
diff --git a/fusl/src/math/sinh.c b/fusl/src/math/sinh.c |
new file mode 100644 |
index 0000000000000000000000000000000000000000..00022c4e6ff6d8250ea238207a12199a1e349dac |
--- /dev/null |
+++ b/fusl/src/math/sinh.c |
@@ -0,0 +1,39 @@ |
+#include "libm.h" |
+ |
+/* sinh(x) = (exp(x) - 1/exp(x))/2 |
+ * = (exp(x)-1 + (exp(x)-1)/exp(x))/2 |
+ * = x + x^3/6 + o(x^5) |
+ */ |
+double sinh(double x) |
+{ |
+ union {double f; uint64_t i;} u = {.f = x}; |
+ uint32_t w; |
+ double t, h, absx; |
+ |
+ h = 0.5; |
+ if (u.i >> 63) |
+ h = -h; |
+ /* |x| */ |
+ u.i &= (uint64_t)-1/2; |
+ absx = u.f; |
+ w = u.i >> 32; |
+ |
+ /* |x| < log(DBL_MAX) */ |
+ if (w < 0x40862e42) { |
+ t = expm1(absx); |
+ if (w < 0x3ff00000) { |
+ if (w < 0x3ff00000 - (26<<20)) |
+ /* note: inexact and underflow are raised by expm1 */ |
+ /* note: this branch avoids spurious underflow */ |
+ return x; |
+ return h*(2*t - t*t/(t+1)); |
+ } |
+ /* note: |x|>log(0x1p26)+eps could be just h*exp(x) */ |
+ return h*(t + t/(t+1)); |
+ } |
+ |
+ /* |x| > log(DBL_MAX) or nan */ |
+ /* note: the result is stored to handle overflow */ |
+ t = 2*h*__expo2(absx); |
+ return t; |
+} |