1  /* @(#)e_atanh.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  /* __ieee754_atanh(x)
      14   * Method :
      15   *    1.Reduced x to positive by atanh(-x) = -atanh(x)
      16   *    2.For x>=0.5
      17   *                  1              2x                          x
      18   *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
      19   *                  2             1 - x                      1 - x
      20   *
      21   *	For x<0.5
      22   *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
      23   *
      24   * Special cases:
      25   *	atanh(x) is NaN if |x| > 1 with signal;
      26   *	atanh(NaN) is that NaN with no signal;
      27   *	atanh(+-1) is +-INF with signal.
      28   *
      29   */
      30  
      31  #include <float.h>
      32  #include <math.h>
      33  #include <math_private.h>
      34  #include <math-underflow.h>
      35  #include <libm-alias-finite.h>
      36  
      37  static const long double one = 1.0L, huge = 1e300L;
      38  
      39  static const long double zero = 0.0L;
      40  
      41  long double
      42  __ieee754_atanhl(long double x)
      43  {
      44  	long double t;
      45  	int64_t hx,ix;
      46  	double xhi;
      47  
      48  	xhi = ldbl_high (x);
      49  	EXTRACT_WORDS64 (hx, xhi);
      50  	ix = hx&0x7fffffffffffffffLL;
      51  	if (ix >= 0x3ff0000000000000LL) { /* |x|>=1 */
      52  	    if (ix > 0x3ff0000000000000LL)
      53  		return (x-x)/(x-x);
      54  	    t = fabsl (x);
      55  	    if (t > one)
      56  		return (x-x)/(x-x);
      57  	    if (t == one)
      58  		return x/zero;
      59  	}
      60  	if(ix<0x3c70000000000000LL&&(huge+x)>zero)	/* x<2**-56 */
      61  	  {
      62  	    math_check_force_underflow (x);
      63  	    return x;
      64  	  }
      65  	x = fabsl (x);
      66  	if(ix<0x3fe0000000000000LL) {		/* x < 0.5 */
      67  	    t = x+x;
      68  	    t = 0.5*__log1pl(t+t*x/(one-x));
      69  	} else
      70  	    t = 0.5*__log1pl((x+x)/(one-x));
      71  	if(hx>=0) return t; else return -t;
      72  }
      73  libm_alias_finite (__ieee754_atanhl, __atanhl)