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