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