1 files changed, 66 insertions, 0 deletions
diff --git a/pl/math/cosh_2u.c b/pl/math/cosh_2u.c
new file mode 100644
index 0000000..5d1df07
--- /dev/null
+++ b/pl/math/cosh_2u.c
@@ -0,0 +1,66 @@
+/*
+ * Double-precision cosh(x) function.
+ *
+ * Copyright (c) 2022-2023, Arm Limited.
+ * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+ */
+
+#include "math_config.h"
+#include "pl_sig.h"
+#include "pl_test.h"
+
+#define AbsMask 0x7fffffffffffffff
+#define SpecialBound                                                           \
+  0x40861da04cbafe44 /* 0x1.61da04cbafe44p+9, above which exp overflows.  */
+
+double
+__exp_dd (double, double);
+
+static double
+specialcase (double x, uint64_t iax)
+{
+  if (iax == 0x7ff0000000000000)
+    return INFINITY;
+  if (iax > 0x7ff0000000000000)
+    return __math_invalid (x);
+  /* exp overflows above SpecialBound. At this magnitude cosh(x) is dominated by
+     exp(x), so we can approximate cosh(x) by (exp(|x|/2)) ^ 2 / 2.  */
+  double t = __exp_dd (asdouble (iax) / 2, 0);
+  return (0.5 * t) * t;
+}
+
+/* Approximation for double-precision cosh(x).
+   cosh(x) = (exp(x) + exp(-x)) / 2.
+   The greatest observed error is in the special region, 1.93 ULP:
+   cosh(0x1.628af341989dap+9) got 0x1.fdf28623ef921p+1021
+			     want 0x1.fdf28623ef923p+1021.
+
+   The greatest observed error in the non-special region is 1.03 ULP:
+   cosh(0x1.502cd8e56ab3bp+0) got 0x1.fe54962842d0ep+0
+			     want 0x1.fe54962842d0fp+0.  */
+double
+cosh (double x)
+{
+  uint64_t ix = asuint64 (x);
+  uint64_t iax = ix & AbsMask;
+
+  /* exp overflows a little bit before cosh, so use special-case handler for the
+     gap, as well as special values.  */
+  if (unlikely (iax >= SpecialBound))
+    return specialcase (x, iax);
+
+  double ax = asdouble (iax);
+  /* Use double-precision exp helper to calculate exp(x), then:
+     cosh(x) = exp(|x|) / 2 + 1 / (exp(|x| * 2).  */
+  double t = __exp_dd (ax, 0);
+  return 0.5 * t + 0.5 / t;
+}
+
+PL_SIG (S, D, 1, cosh, -10.0, 10.0)
+PL_TEST_ULP (cosh, 1.43)
+PL_TEST_INTERVAL (cosh, 0, 0x1.61da04cbafe44p+9, 100000)
+PL_TEST_INTERVAL (cosh, -0, -0x1.61da04cbafe44p+9, 100000)
+PL_TEST_INTERVAL (cosh, 0x1.61da04cbafe44p+9, 0x1p10, 1000)
+PL_TEST_INTERVAL (cosh, -0x1.61da04cbafe44p+9, -0x1p10, 1000)
+PL_TEST_INTERVAL (cosh, 0x1p10, inf, 100)
+PL_TEST_INTERVAL (cosh, -0x1p10, -inf, 100)