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_inlines_amd.h"
 31  #include "libm_special.h"
 32  
 33  double FN_PROTOTYPE_REF(tanh)(double x)
 34  {
 35    /*
 36      The definition of tanh(x) is sinh(x)/cosh(x), which is also equivalent
 37      to the following three formulae:
 38      1.  (exp(x) - exp(-x))/(exp(x) + exp(-x))
 39      2.  (1 - (2/(exp(2*x) + 1 )))
 40      3.  (exp(2*x) - 1)/(exp(2*x) + 1)
 41      but computationally, some formulae are better on some ranges.
 42    */
 43    static const double
 44      thirtytwo_by_log2 = 4.61662413084468283841e+01, /* 0x40471547652b82fe */
 45      log2_by_32_lead = 2.16608493356034159660e-02, /* 0x3f962e42fe000000 */
 46      log2_by_32_tail = 5.68948749532545630390e-11, /* 0x3dcf473de6af278e */
 47      large_threshold = 20.0; /* 0x4034000000000000 */
 48  
 49    unsigned long long ux, aux, xneg;
 50    double y, z, p, z1, z2;
 51    int m;
 52  
 53    /* Special cases */
 54  
 55    GET_BITS_DP64(x, ux);
 56    aux = ux & ~SIGNBIT_DP64;
 57    if (aux < 0x3e30000000000000) /* |x| small enough that tanh(x) = x */
 58      {      
 59          if (aux == 0x0000000000000000)
 60          {
 61              /* x is +/-zero. Return the same zero. */
 62              return x;
 63          }
 64          else
 65          {
 66  #ifdef WINDOWS
 67              return x;
 68  #else
 69              return __amd_handle_error("tanh", __amd_tanh, ux, _UNDERFLOW, AMD_F_INEXACT|AMD_F_UNDERFLOW, ERANGE, x, 0.0, 1);
 70  #endif
 71          }
 72      }
 73    else if  (aux > 0x7ff0000000000000) /* |x| is NaN */
 74      {
 75  #ifdef WINDOWS
 76          return __amd_handle_error("tanh", __amd_tanh, ux|QNANBITPATT_DP64, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0, 1);
 77  #else
 78          return x+x;
 79  #endif
 80      }
 81  
 82    xneg = (aux != ux);
 83  
 84    y = x;
 85    if (xneg) y = -x;
 86  
 87    if (y > large_threshold)
 88      {
 89        /* If x is large then exp(-x) is negligible and
 90           formula 1 reduces to plus or minus 1.0 */
 91        z = 1.0;
 92      }
 93    else if (y <= 1.0)
 94      {
 95        double y2;
 96        y2 = y*y;
 97        if (y < 0.9)
 98          {
 99            /* Use a [3,3] Remez approximation on [0,0.9]. */
100            z = y + y*y2*
101              (-0.274030424656179760118928e0 +
102               (-0.176016349003044679402273e-1 +
103                (-0.200047621071909498730453e-3 -
104                 0.142077926378834722618091e-7*y2)*y2)*y2)/
105              (0.822091273968539282568011e0 +
106               (0.381641414288328849317962e0 +
107                (0.201562166026937652780575e-1 +
108                 0.2091140262529164482568557e-3*y2)*y2)*y2);
109          }
110        else
111          {
112            /* Use a [3,3] Remez approximation on [0.9,1]. */
113            z = y + y*y2*
114              (-0.227793870659088295252442e0 +
115               (-0.146173047288731678404066e-1 +
116                (-0.165597043903549960486816e-3 -
117                 0.115475878996143396378318e-7*y2)*y2)*y2)/
118              (0.683381611977295894959554e0 +
119               (0.317204558977294374244770e0 +
120                (0.167358775461896562588695e-1 +
121                 0.173076050126225961768710e-3*y2)*y2)*y2);
122          }
123      }
124    else
125      {
126        /* Compute p = exp(2*y) + 1. The code is basically inlined
127           from exp_amd. */
128  
129        splitexp(2*y, 1.0, thirtytwo_by_log2, log2_by_32_lead,
130             log2_by_32_tail, &m, &z1, &z2);
131        p = scaleDouble_2(z1 + z2, m) + 1.0;
132  
133        /* Now reconstruct tanh from p. */
134        z = (1.0 - 2.0/p);
135      }
136  
137    if (xneg) z = - z;
138    return z;
139  }
140  
141  
142