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-barriers.h>
      36  #include <math_private.h>
      37  #include <math-underflow.h>
      38  #include <libm-alias-finite.h>
      39  
      40  static const long double one = 1.0, huge = 1e4900L;
      41  
      42  static const long double zero = 0.0;
      43  
      44  long double
      45  __ieee754_atanhl(long double x)
      46  {
      47  	long double t;
      48  	int32_t ix;
      49  	uint32_t se,i0,i1;
      50  	GET_LDOUBLE_WORDS(se,i0,i1,x);
      51  	ix = se&0x7fff;
      52  	if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff)
      53  	  /* |x|>1 */
      54  	    return (x-x)/(x-x);
      55  	if(ix==0x3fff)
      56  	    return x/zero;
      57  	if(ix<0x3fdf) {
      58  	    math_force_eval(huge+x);
      59  	    math_check_force_underflow (x);
      60  	    return x;	/* x<2**-32 */
      61  	}
      62  	SET_LDOUBLE_EXP(x,ix);
      63  	if(ix<0x3ffe) {		/* x < 0.5 */
      64  	    t = x+x;
      65  	    t = 0.5*__log1pl(t+t*x/(one-x));
      66  	} else
      67  	    t = 0.5*__log1pl((x+x)/(one-x));
      68  	if(se<=0x7fff) return t; else return -t;
      69  }
      70  libm_alias_finite (__ieee754_atanhl, __atanhl)