| Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame^] | 1 | /* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_asinhl.c */ |
| 2 | /* |
| 3 | * ==================================================== |
| 4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 5 | * |
| 6 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 7 | * Permission to use, copy, modify, and distribute this |
| 8 | * software is freely granted, provided that this notice |
| 9 | * is preserved. |
| 10 | * ==================================================== |
| 11 | */ |
| 12 | /* asinhl(x) |
| 13 | * Method : |
| 14 | * Based on |
| 15 | * asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ] |
| 16 | * we have |
| 17 | * asinhl(x) := x if 1+x*x=1, |
| 18 | * := signl(x)*(logl(x)+ln2)) for large |x|, else |
| 19 | * := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else |
| 20 | * := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2))) |
| 21 | */ |
| 22 | |
| 23 | #include "libm.h" |
| 24 | |
| 25 | #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| 26 | long double asinhl(long double x) |
| 27 | { |
| 28 | return asinh(x); |
| 29 | } |
| 30 | #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
| 31 | static const long double |
| 32 | one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */ |
| 33 | ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ |
| 34 | huge = 1.000000000000000000e+4900L; |
| 35 | |
| 36 | long double asinhl(long double x) |
| 37 | { |
| 38 | long double t,w; |
| 39 | int32_t hx,ix; |
| 40 | |
| 41 | GET_LDOUBLE_EXP(hx, x); |
| 42 | ix = hx & 0x7fff; |
| 43 | if (ix == 0x7fff) |
| 44 | return x + x; /* x is inf or NaN */ |
| 45 | if (ix < 0x3fde) { /* |x| < 2**-34 */ |
| 46 | /* return x, raise inexact if x != 0 */ |
| 47 | if (huge+x > one) |
| 48 | return x; |
| 49 | } |
| 50 | if (ix > 0x4020) { /* |x| > 2**34 */ |
| 51 | w = logl(fabsl(x)) + ln2; |
| 52 | } else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */ |
| 53 | t = fabsl(x); |
| 54 | w = logl(2.0*t + one/(sqrtl(x*x + one) + t)); |
| 55 | } else { /* 2.0 > |x| > 2**-28 */ |
| 56 | t = x*x; |
| 57 | w =log1pl(fabsl(x) + t/(one + sqrtl(one + t))); |
| 58 | } |
| 59 | if (hx & 0x8000) |
| 60 | return -w; |
| 61 | return w; |
| 62 | } |
| 63 | #endif |