1  /*
  2   * Copyright (C) 2008-2020 Advanced Micro Devices, Inc. All rights reserved.
  3   *
  4   * Redistribution and use in source and binary forms, with or without modification,
  5   * are permitted provided that the following conditions are met:
  6   * 1. Redistributions of source code must retain the above copyright notice,
  7   *    this list of conditions and the following disclaimer.
  8   * 2. Redistributions in binary form must reproduce the above copyright notice,
  9   *    this list of conditions and the following disclaimer in the documentation
 10   *    and/or other materials provided with the distribution.
 11   * 3. Neither the name of the copyright holder nor the names of its contributors
 12   *    may be used to endorse or promote products derived from this software without
 13   *    specific prior written permission.
 14   *
 15   * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 16   * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 17   * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 18   * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
 19   * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 20   * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 21   * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 22   * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 23   * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 24   * POSSIBILITY OF SUCH DAMAGE.
 25   *
 26   */
 27  
 28  #include "libm_amd.h"
 29  #include "libm_util_amd.h"
 30  #include "libm_special.h"
 31  
 32  #undef _FUNCNAME
 33  #define _FUNCNAME "atanhf"
 34  float FN_PROTOTYPE_REF(atanhf)(float x)
 35  {
 36  
 37    double dx;
 38    unsigned int ux, ax;
 39    double r, t, poly;
 40  
 41    GET_BITS_SP32(x, ux);
 42    ax = ux & ~SIGNBIT_SP32;
 43  
 44    if ((ux & EXPBITS_SP32) == EXPBITS_SP32)
 45      {
 46        /* x is either NaN or infinity */
 47        if (ux & MANTBITS_SP32)
 48          {
 49            /* x is NaN */
 50  #ifdef WINDOWS
 51  	return __amd_handle_errorf(_FUNCNAME,__amd_tanh, ux|0x00400000, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0F,1);
 52  #else
 53  	if (ux & QNAN_MASK_32)
 54  	return __amd_handle_errorf(_FUNCNAME,__amd_tanh, ux|0x00400000, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0F,1);
 55  	else
 56  	return __amd_handle_errorf(_FUNCNAME,__amd_tanh, ux|0x00400000, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0F,1);
 57  #endif
 58          }
 59        else
 60          {
 61            /* x is infinity; return a NaN */
 62  			return __amd_handle_errorf(_FUNCNAME,__amd_tanh, INDEFBITPATT_SP32, _DOMAIN,
 63                                AMD_F_INVALID, EDOM, x, 0.0F,1);
 64  
 65          }
 66      }
 67    else if (ax >= 0x3f800000)
 68      {
 69        if (ax > 0x3f800000)
 70          {
 71            /* abs(x) > 1.0; return NaN */
 72            return __amd_handle_errorf(_FUNCNAME,__amd_tanh, INDEFBITPATT_SP32, _DOMAIN,
 73                                 AMD_F_INVALID, EDOM, x, 0.0F,1);
 74          }
 75        else if (ux == 0x3f800000)
 76          {
 77            /* x = +1.0; return infinity with the same sign as x
 78               and set the divbyzero status flag */
 79            return __amd_handle_errorf(_FUNCNAME,__amd_tanh, PINFBITPATT_SP32, _DOMAIN,
 80                                 AMD_F_DIVBYZERO, EDOM, x, 0.0F,1);
 81          }
 82        else
 83          {
 84            /* x = -1.0; return infinity with the same sign as x */
 85            return __amd_handle_errorf(_FUNCNAME,__amd_tanh, NINFBITPATT_SP32, _DOMAIN,
 86                                 AMD_F_DIVBYZERO, EDOM, x, 0.0F,1);
 87          }
 88      }
 89  
 90    if (ax < 0x39000000)
 91      {
 92        if (ax == 0x00000000)
 93          {
 94            /* x is +/-zero. Return the same zero. */
 95            return x;
 96          }
 97        else
 98          {
 99            /* Arguments smaller than 2^(-13) in magnitude are
100               approximated by atanhf(x) = x, raising inexact flag. */
101  #ifdef WINDOWS
102            return x; // return valf_with_flags(x, AMD_F_INEXACT);
103  #else
104  	return __amd_handle_errorf(_FUNCNAME,__amd_tanh, ux, _UNDERFLOW, AMD_F_UNDERFLOW|AMD_F_INEXACT, ERANGE, x, 0.0F,1);
105  #endif
106          }
107      }
108    else
109      {
110        dx = x;
111        if (ax < 0x3f000000)
112          {
113            /* Arguments up to 0.5 in magnitude are
114               approximated by a [2,2] minimax polynomial */
115            t = dx*dx;
116            poly =
117              (0.39453629046e0 +
118             (-0.28120347286e0 +
119               0.92834212715e-2 * t) * t) /
120              (0.11836088638e1 + 
121             (-0.15537744551e1 +
122               0.45281890445e0 * t) * t);
123            return (float)(dx + dx*t*poly);
124          }
125        else
126          {
127            /* abs(x) >= 0.5 */
128            /* Note that
129                 atanhf(x) = 0.5 * ln((1+x)/(1-x))
130               (see Abramowitz and Stegun 4.6.22).
131               For greater accuracy we use the variant formula
132               atanhf(x) = log(1 + 2x/(1-x)) = log1p(2x/(1-x)).
133            */
134            if (ux & SIGNBIT_SP32)
135              {
136                /* Argument x is negative */
137                r = (-2.0 * dx) / (1.0 + dx);
138                r = 0.5 * FN_PROTOTYPE(log1p)(r);
139                return (float)-r;
140              }
141            else
142              {
143                r = (2.0 * dx) / (1.0 - dx);
144                r = 0.5 * FN_PROTOTYPE(log1p)(r);
145                return (float)r;
146              }
147          }
148      }
149  }
150  
151