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.c58
1 files changed, 12 insertions, 46 deletions
diff --git a/src/math/acoshl.c b/src/math/acoshl.c
index a4024516..472c71cb 100644
--- a/src/math/acoshl.c
+++ b/src/math/acoshl.c
@@ -1,28 +1,3 @@
-/* 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
@@ -31,29 +6,20 @@ long double acoshl(long double x)
return acosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double
-ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
-
+/* acosh(x) = log(x + sqrt(x*x-1)) */
long double acoshl(long double x)
{
- long double t;
- uint32_t se,i0,i1;
+ union {
+ long double f;
+ struct{uint64_t m; int16_t se; uint16_t pad;} i;
+ } u = {.f = x};
- 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 - 1.0/(x + sqrtl(t - 1.0)));
- }
- /* 1 < x <= 2 */
- t = x - 1.0;
- return log1pl(t + sqrtl(2.0*t + t*t));
+ if (u.i.se < 0x3fff + 1)
+ /* x < 2, invalid if x < 1 or nan */
+ return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1)));
+ if (u.i.se < 0x3fff + 32)
+ /* x < 0x1p32 */
+ return logl(2*x - 1/(x+sqrtl(x*x-1)));
+ return logl(x) + 0.693147180559945309417232121458176568L;
}
#endif