1/* s_atanhl.c -- __float128 version of s_atan.c. 2 * Conversion to __float128 by Ulrich Drepper, 3 * Cygnus Support, drepper@cygnus.com. 4 */ 5 6/* 7 * ==================================================== 8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 9 * 10 * Developed at SunPro, a Sun Microsystems, Inc. business. 11 * Permission to use, copy, modify, and distribute this 12 * software is freely granted, provided that this notice 13 * is preserved. 14 * ==================================================== 15 */ 16 17/* atanhq(x) 18 * Method : 19 * 1.Reduced x to positive by atanh(-x) = -atanh(x) 20 * 2.For x>=0.5 21 * 1 2x x 22 * atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 23 * 2 1 - x 1 - x 24 * 25 * For x<0.5 26 * atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x)) 27 * 28 * Special cases: 29 * atanhl(x) is NaN if |x| > 1 with signal; 30 * atanhl(NaN) is that NaN with no signal; 31 * atanhl(+-1) is +-INF with signal. 32 * 33 */ 34 35#include "quadmath-imp.h" 36 37static const __float128 one = 1.0Q, huge = 1e4900Q; 38static const __float128 zero = 0.0Q; 39 40__float128 41atanhq (__float128 x) 42{ 43 __float128 t; 44 uint32_t jx, ix; 45 ieee854_float128 u; 46 47 u.value = x; 48 jx = u.words32.w0; 49 ix = jx & 0x7fffffff; 50 u.words32.w0 = ix; 51 if (ix >= 0x3fff0000) /* |x| >= 1.0 or infinity or NaN */ 52 { 53 if (u.value == one) 54 return x/zero; 55 else 56 return (x-x)/(x-x); 57 } 58 if(ix<0x3fc60000 && (huge+x)>zero) return x; /* x < 2^-57 */ 59 60 if(ix<0x3ffe0000) { /* x < 0.5 */ 61 t = u.value+u.value; 62 t = 0.5*log1pq(t+t*u.value/(one-u.value)); 63 } else 64 t = 0.5*log1pq((u.value+u.value)/(one-u.value)); 65 if(jx & 0x80000000) return -t; else return t; 66} 67