1*49393c00Smartynas /* @(#)e_atanh.c 5.1 93/09/24 */ 2*49393c00Smartynas /* 3*49393c00Smartynas * ==================================================== 4*49393c00Smartynas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5*49393c00Smartynas * 6*49393c00Smartynas * Developed at SunPro, a Sun Microsystems, Inc. business. 7*49393c00Smartynas * Permission to use, copy, modify, and distribute this 8*49393c00Smartynas * software is freely granted, provided that this notice 9*49393c00Smartynas * is preserved. 10*49393c00Smartynas * ==================================================== 11*49393c00Smartynas */ 12*49393c00Smartynas 13*49393c00Smartynas /* atanhl(x) 14*49393c00Smartynas * Method : 15*49393c00Smartynas * 1.Reduced x to positive by atanh(-x) = -atanh(x) 16*49393c00Smartynas * 2.For x>=0.5 17*49393c00Smartynas * 1 2x x 18*49393c00Smartynas * atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 19*49393c00Smartynas * 2 1 - x 1 - x 20*49393c00Smartynas * 21*49393c00Smartynas * For x<0.5 22*49393c00Smartynas * atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x)) 23*49393c00Smartynas * 24*49393c00Smartynas * Special cases: 25*49393c00Smartynas * atanhl(x) is NaN if |x| > 1 with signal; 26*49393c00Smartynas * atanhl(NaN) is that NaN with no signal; 27*49393c00Smartynas * atanhl(+-1) is +-INF with signal. 28*49393c00Smartynas * 29*49393c00Smartynas */ 30*49393c00Smartynas 31*49393c00Smartynas #include <math.h> 32*49393c00Smartynas 33*49393c00Smartynas #include "math_private.h" 34*49393c00Smartynas 35*49393c00Smartynas static const long double one = 1.0L, huge = 1e4900L; 36*49393c00Smartynas 37*49393c00Smartynas static const long double zero = 0.0L; 38*49393c00Smartynas 39*49393c00Smartynas long double atanhl(long double x)40*49393c00Smartynasatanhl(long double x) 41*49393c00Smartynas { 42*49393c00Smartynas long double t; 43*49393c00Smartynas u_int32_t jx, ix; 44*49393c00Smartynas ieee_quad_shape_type u; 45*49393c00Smartynas 46*49393c00Smartynas u.value = x; 47*49393c00Smartynas jx = u.parts32.mswhi; 48*49393c00Smartynas ix = jx & 0x7fffffff; 49*49393c00Smartynas u.parts32.mswhi = ix; 50*49393c00Smartynas if (ix >= 0x3fff0000) /* |x| >= 1.0 or infinity or NaN */ 51*49393c00Smartynas { 52*49393c00Smartynas if (u.value == one) 53*49393c00Smartynas return x/zero; 54*49393c00Smartynas else 55*49393c00Smartynas return (x-x)/(x-x); 56*49393c00Smartynas } 57*49393c00Smartynas if(ix<0x3fc60000 && (huge+x)>zero) return x; /* x < 2^-57 */ 58*49393c00Smartynas 59*49393c00Smartynas if(ix<0x3ffe0000) { /* x < 0.5 */ 60*49393c00Smartynas t = u.value+u.value; 61*49393c00Smartynas t = 0.5*log1pl(t+t*u.value/(one-u.value)); 62*49393c00Smartynas } else 63*49393c00Smartynas t = 0.5*log1pl((u.value+u.value)/(one-u.value)); 64*49393c00Smartynas if(jx & 0x80000000) return -t; else return t; 65*49393c00Smartynas } 66