1  /* Copyright (C) 2011-2023 Free Software Foundation, Inc.
       2     This file is part of the GNU C Library.
       3  
       4     The GNU C Library is free software; you can redistribute it and/or
       5     modify it under the terms of the GNU Lesser General Public
       6     License as published by the Free Software Foundation; either
       7     version 2.1 of the License, or (at your option) any later version.
       8  
       9     The GNU C Library is distributed in the hope that it will be useful,
      10     but WITHOUT ANY WARRANTY; without even the implied warranty of
      11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12     Lesser General Public License for more details.
      13  
      14     You should have received a copy of the GNU Lesser General Public
      15     License along with the GNU C Library; if not, see
      16     <https://www.gnu.org/licenses/>.  */
      17  
      18  
      19  /* __ieee754_atanh(x)
      20     Method :
      21        1.Reduced x to positive by atanh(-x) = -atanh(x)
      22        2.For x>=0.5
      23  		    1              2x                          x
      24  	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
      25  		    2             1 - x                      1 - x
      26  
      27  	For x<0.5
      28  	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
      29  
      30     Special cases:
      31  	atanh(x) is NaN if |x| > 1 with signal;
      32  	atanh(NaN) is that NaN with no signal;
      33  	atanh(+-1) is +-INF with signal.
      34  
      35   */
      36  
      37  #include <float.h>
      38  #include <inttypes.h>
      39  #include <math.h>
      40  #include <math-barriers.h>
      41  #include <math_private.h>
      42  #include <math-underflow.h>
      43  #include <libm-alias-finite.h>
      44  
      45  static const float huge = 1e30;
      46  
      47  float
      48  __ieee754_atanhf (float x)
      49  {
      50    float xa = fabsf (x);
      51    float t;
      52    if (isless (xa, 0.5f))
      53      {
      54        if (__glibc_unlikely (xa < 0x1.0p-28f))
      55  	{
      56  	  math_force_eval (huge + x);
      57  	  math_check_force_underflow (x);
      58  	  return x;
      59  	}
      60  
      61        t = xa + xa;
      62        t = 0.5f * __log1pf (t + t * xa / (1.0f - xa));
      63      }
      64    else if (__glibc_likely (isless (xa, 1.0f)))
      65      t = 0.5f * __log1pf ((xa + xa) / (1.0f - xa));
      66    else
      67      {
      68        if (isgreater (xa, 1.0f))
      69  	return (x - x) / (x - x);
      70  
      71        return x / 0.0f;
      72      }
      73  
      74    return copysignf (t, x);
      75  }
      76  libm_alias_finite (__ieee754_atanhf, __atanhf)