summaryrefslogtreecommitdiff
path: root/src/math/acoshl.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/math/acoshl.c')
-rw-r--r--src/math/acoshl.c60
1 files changed, 60 insertions, 0 deletions
diff --git a/src/math/acoshl.c b/src/math/acoshl.c
new file mode 100644
index 00000000..d8310a73
--- /dev/null
+++ b/src/math/acoshl.c
@@ -0,0 +1,60 @@
+/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_acoshl.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.
+ * ====================================================
+ */
+/* acoshl(x)
+ * Method :
+ * Based on
+ * acoshl(x) = logl [ x + sqrtl(x*x-1) ]
+ * we have
+ * acoshl(x) := logl(x)+ln2, if x is large; else
+ * acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
+ * acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ * acoshl(x) is NaN with signal if x<1.
+ * acoshl(NaN) is NaN without signal.
+ */
+
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+long double acoshl(long double x)
+{
+ return acosh(x);
+}
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+static const long double
+one = 1.0,
+ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
+
+long double acoshl(long double x)
+{
+ long double t;
+ uint32_t se,i0,i1;
+
+ GET_LDOUBLE_WORDS(se, i0, i1, x);
+ if (se < 0x3fff || se & 0x8000) { /* x < 1 */
+ return (x-x)/(x-x);
+ } else if (se >= 0x401d) { /* x > 2**30 */
+ if (se >= 0x7fff) /* x is inf or NaN */
+ return x+x;
+ return logl(x) + ln2; /* acoshl(huge) = logl(2x) */
+ } else if (((se-0x3fff)|i0|i1) == 0) {
+ return 0.0; /* acosh(1) = 0 */
+ } else if (se > 0x4000) { /* x > 2 */
+ t = x*x;
+ return logl(2.0*x - one/(x + sqrtl(t - one)));
+ }
+ /* 1 < x <= 2 */
+ t = x - one;
+ return log1pl(t + sqrtl(2.0*t + t*t));
+}
+#endif