1 /* @(#)e_acosh.c 5.1 93/09/24 */ 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 13 /* acoshl(x) 14 * Method : 15 * Based on 16 * acoshl(x) = logl [ x + sqrtl(x*x-1) ] 17 * we have 18 * acoshl(x) := logl(x)+ln2, if x is large; else 19 * acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else 20 * acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1. 21 * 22 * Special cases: 23 * acoshl(x) is NaN with signal if x<1. 24 * acoshl(NaN) is NaN without signal. 25 */ 26 27 #include <openlibm_math.h> 28 29 #include "math_private.h" 30 31 static const long double 32 one = 1.0, 33 ln2 = 0.6931471805599453094172321214581766L; 34 35 long double 36 acoshl(long double x) 37 { 38 long double t; 39 u_int64_t lx; 40 int64_t hx; 41 GET_LDOUBLE_WORDS64(hx,lx,x); 42 if(hx<0x3fff000000000000LL) { /* x < 1 */ 43 return (x-x)/(x-x); 44 } else if(hx >=0x4035000000000000LL) { /* x > 2**54 */ 45 if(hx >=0x7fff000000000000LL) { /* x is inf of NaN */ 46 return x+x; 47 } else 48 return logl(x)+ln2; /* acoshl(huge)=logl(2x) */ 49 } else if(((hx-0x3fff000000000000LL)|lx)==0) { 50 return 0.0L; /* acosh(1) = 0 */ 51 } else if (hx > 0x4000000000000000LL) { /* 2**28 > x > 2 */ 52 t=x*x; 53 return logl(2.0L*x-one/(x+sqrtl(t-one))); 54 } else { /* 1<x<2 */ 55 t = x-one; 56 return log1pl(t+sqrtl(2.0L*t+t*t)); 57 } 58 } 59