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 "asinhf"
 34  float FN_PROTOTYPE_REF(asinhf)(float x)
 35  {
 36  
 37    double dx;
 38    unsigned int ux, ax, xneg;
 39    double absx, r, rarg, t, poly;
 40  
 41    static const unsigned int
 42      rteps = 0x39800000,    /* sqrt(eps) = 2.44140625000000000000e-04 */
 43      recrteps = 0x46000000; /* 1/rteps = 4.09600000000000000000e+03 */
 44  
 45    static const double
 46      log2 = 6.93147180559945286227e-01;  /* 0x3fe62e42fefa39ef */
 47  
 48    GET_BITS_SP32(x, ux);
 49    ax = ux & ~SIGNBIT_SP32;
 50    xneg = ux & SIGNBIT_SP32;
 51  
 52    if ((ux & EXPBITS_SP32) == EXPBITS_SP32)
 53      {
 54        /* x is either NaN or infinity */
 55        if (ux & MANTBITS_SP32)
 56          {
 57            /* x is NaN */
 58  #ifdef WINDOWS
 59            return __amd_handle_errorf(_FUNCNAME,__amd_asinh, ux|0x00400000, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0F,1);
 60  #else
 61  	  if (ux & QNAN_MASK_32)
 62            return __amd_handle_errorf(_FUNCNAME,__amd_asinh, ux|0x00400000, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0F,1);
 63  	  else
 64            return __amd_handle_errorf(_FUNCNAME,__amd_asinh, ux|0x00400000, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0F,1);
 65  #endif
 66          }
 67        else
 68          {
 69            /* x is infinity. Return the same infinity. */
 70  			return x;
 71          }
 72      }
 73    else if (ax < rteps) /* abs(x) < sqrt(epsilon) */
 74      {
 75        if (ax == 0x00000000)
 76          {
 77            /* x is +/-zero. Return the same zero. */
 78            return x;
 79          }
 80        else
 81          {
 82            /* Tiny arguments approximated by asinhf(x) = x
 83               - avoid slow operations on denormalized numbers */
 84  #ifdef WINDOWS
 85  			return x; //return valf_with_flags(x,AMD_F_INEXACT);
 86  #else
 87            return __amd_handle_errorf(_FUNCNAME,__amd_asinh, ux, _UNDERFLOW, AMD_F_UNDERFLOW|AMD_F_INEXACT, ERANGE, x, 0.0F,1);
 88  #endif
 89          }
 90      }
 91  
 92    dx = x;
 93    if (xneg)
 94      absx = -dx;
 95    else
 96      absx = dx;
 97  
 98    if (ax <= 0x40800000) /* abs(x) <= 4.0 */
 99      {
100        /* Arguments less than 4.0 in magnitude are
101           approximated by [4,4] minimax polynomials
102        */
103        t = dx*dx;
104        if (ax <= 0x40000000) /* abs(x) <= 2 */
105          poly =
106            (-0.1152965835871758072e-1 +
107            (-0.1480204186473758321e-1 +
108            (-0.5063201055468483248e-2 +
109            (-0.4162727710583425360e-3 -
110              0.1177198915954942694e-5 * t) * t) * t) * t) /
111             (0.6917795026025976739e-1 +
112             (0.1199423176003939087e+0 +
113             (0.6582362487198468066e-1 +
114             (0.1260024978680227945e-1 +
115              0.6284381367285534560e-3 * t) * t) * t) * t);
116        else
117          poly =
118             (-0.185462290695578589e-2 +
119             (-0.113672533502734019e-2 +
120             (-0.142208387300570402e-3 +
121             (-0.339546014993079977e-5 -
122               0.151054665394480990e-8 * t) * t) * t) * t) /
123              (0.111486158580024771e-1 +
124              (0.117782437980439561e-1 +
125              (0.325903773532674833e-2 +
126              (0.255902049924065424e-3 +
127               0.434150786948890837e-5 * t) * t) * t) * t);
128        return (float)(dx + dx*t*poly);
129      }
130    else
131      {
132        /* abs(x) > 4.0 */
133        if (ax > recrteps)
134          {
135            /* Arguments greater than 1/sqrt(epsilon) in magnitude are
136               approximated by asinhf(x) = ln(2) + ln(abs(x)), with sign of x */
137            r = FN_PROTOTYPE(log)(absx) + log2;
138          }
139        else
140          {
141            rarg = absx*absx+1.0;
142            /* Arguments such that 4.0 <= abs(x) <= 1/sqrt(epsilon) are
143               approximated by
144                 asinhf(x) = ln(abs(x) + sqrt(x*x+1))
145               with the sign of x (see Abramowitz and Stegun 4.6.20) */
146            /* Use assembly instruction to compute r = sqrt(rarg); */
147            ASMSQRT(rarg,r);
148            r += absx;
149            r = FN_PROTOTYPE(log)(r);
150          }
151        if (xneg)
152          return (float)(-r);
153        else
154          return (float)r;
155      }
156  }
157  
158