1*181254a7Smrg /* e_acoshl.c -- long double version of e_acosh.c.
2*181254a7Smrg * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
3*181254a7Smrg */
4*181254a7Smrg
5*181254a7Smrg /*
6*181254a7Smrg * ====================================================
7*181254a7Smrg * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8*181254a7Smrg *
9*181254a7Smrg * Developed at SunPro, a Sun Microsystems, Inc. business.
10*181254a7Smrg * Permission to use, copy, modify, and distribute this
11*181254a7Smrg * software is freely granted, provided that this notice
12*181254a7Smrg * is preserved.
13*181254a7Smrg * ====================================================
14*181254a7Smrg */
15*181254a7Smrg
16*181254a7Smrg /* acoshq(x)
17*181254a7Smrg * Method :
18*181254a7Smrg * Based on
19*181254a7Smrg * acoshl(x) = logq [ x + sqrtq(x*x-1) ]
20*181254a7Smrg * we have
21*181254a7Smrg * acoshl(x) := logq(x)+ln2, if x is large; else
22*181254a7Smrg * acoshl(x) := logq(2x-1/(sqrtq(x*x-1)+x)) if x>2; else
23*181254a7Smrg * acoshl(x) := log1pq(t+sqrtq(2.0*t+t*t)); where t=x-1.
24*181254a7Smrg *
25*181254a7Smrg * Special cases:
26*181254a7Smrg * acoshl(x) is NaN with signal if x<1.
27*181254a7Smrg * acoshl(NaN) is NaN without signal.
28*181254a7Smrg */
29*181254a7Smrg
30*181254a7Smrg #include "quadmath-imp.h"
31*181254a7Smrg
32*181254a7Smrg static const __float128
33*181254a7Smrg one = 1.0,
34*181254a7Smrg ln2 = 0.6931471805599453094172321214581766Q;
35*181254a7Smrg
36*181254a7Smrg __float128
acoshq(__float128 x)37*181254a7Smrg acoshq(__float128 x)
38*181254a7Smrg {
39*181254a7Smrg __float128 t;
40*181254a7Smrg uint64_t lx;
41*181254a7Smrg int64_t hx;
42*181254a7Smrg GET_FLT128_WORDS64(hx,lx,x);
43*181254a7Smrg if(hx<0x3fff000000000000LL) { /* x < 1 */
44*181254a7Smrg return (x-x)/(x-x);
45*181254a7Smrg } else if(hx >=0x4035000000000000LL) { /* x > 2**54 */
46*181254a7Smrg if(hx >=0x7fff000000000000LL) { /* x is inf of NaN */
47*181254a7Smrg return x+x;
48*181254a7Smrg } else
49*181254a7Smrg return logq(x)+ln2; /* acoshl(huge)=logq(2x) */
50*181254a7Smrg } else if(((hx-0x3fff000000000000LL)|lx)==0) {
51*181254a7Smrg return 0; /* acosh(1) = 0 */
52*181254a7Smrg } else if (hx > 0x4000000000000000LL) { /* 2**28 > x > 2 */
53*181254a7Smrg t=x*x;
54*181254a7Smrg return logq(2*x-one/(x+sqrtq(t-one)));
55*181254a7Smrg } else { /* 1<x<2 */
56*181254a7Smrg t = x-one;
57*181254a7Smrg return log1pq(t+sqrtq(2*t+t*t));
58*181254a7Smrg }
59*181254a7Smrg }
60