Index: fusl/src/math/cosh.c |
diff --git a/fusl/src/math/cosh.c b/fusl/src/math/cosh.c |
new file mode 100644 |
index 0000000000000000000000000000000000000000..100f8231d8829604b9e293491a9abb1774b4b29a |
--- /dev/null |
+++ b/fusl/src/math/cosh.c |
@@ -0,0 +1,40 @@ |
+#include "libm.h" |
+ |
+/* cosh(x) = (exp(x) + 1/exp(x))/2 |
+ * = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x) |
+ * = 1 + x*x/2 + o(x^4) |
+ */ |
+double cosh(double x) |
+{ |
+ union {double f; uint64_t i;} u = {.f = x}; |
+ uint32_t w; |
+ double t; |
+ |
+ /* |x| */ |
+ u.i &= (uint64_t)-1/2; |
+ x = u.f; |
+ w = u.i >> 32; |
+ |
+ /* |x| < log(2) */ |
+ if (w < 0x3fe62e42) { |
+ if (w < 0x3ff00000 - (26<<20)) { |
+ /* raise inexact if x!=0 */ |
+ FORCE_EVAL(x + 0x1p120f); |
+ return 1; |
+ } |
+ t = expm1(x); |
+ return 1 + t*t/(2*(1+t)); |
+ } |
+ |
+ /* |x| < log(DBL_MAX) */ |
+ if (w < 0x40862e42) { |
+ t = exp(x); |
+ /* note: if x>log(0x1p26) then the 1/t is not needed */ |
+ return 0.5*(t + 1/t); |
+ } |
+ |
+ /* |x| > log(DBL_MAX) or nan */ |
+ /* note: the result is stored to handle overflow */ |
+ t = __expo2(x); |
+ return t; |
+} |