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  
 33  #undef _FUNCNAME
 34  #define _FUNCNAME "acoshf"
 35  float FN_PROTOTYPE_REF(acoshf)(float x)
 36  {
 37  
 38    unsigned int ux;
 39    double dx, r, rarg, t;
 40  
 41    static const unsigned int
 42      recrteps = 0x46000000; /* 1/sqrt(eps) = 4.09600000000000000000e+03 */
 43  
 44    static const double
 45      log2 = 6.93147180559945286227e-01;  /* 0x3fe62e42fefa39ef */
 46  
 47    GET_BITS_SP32(x, ux);
 48  
 49    if ((ux & EXPBITS_SP32) == EXPBITS_SP32)
 50      {
 51        /* x is either NaN or infinity */
 52        if (ux & MANTBITS_SP32)
 53          {
 54            /* x is NaN */
 55  #ifdef WINDOWS
 56            return __amd_handle_errorf(_FUNCNAME,__amd_acosh, ux|0x00400000, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0F,1);
 57  #else
 58  		if(ux & QNAN_MASK_32)
 59            return __amd_handle_errorf(_FUNCNAME,__amd_acosh, ux|0x00400000, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0F,1);
 60  		else
 61            return __amd_handle_errorf(_FUNCNAME,__amd_acosh, ux|0x00400000, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0F,1);
 62  #endif
 63         }
 64        else
 65          {
 66            /* x is infinity */
 67            if (ux & SIGNBIT_SP32)
 68              /* x is negative infinity. Return a NaN. */
 69              return __amd_handle_errorf(_FUNCNAME,__amd_acosh, INDEFBITPATT_SP32, _DOMAIN,
 70                                   AMD_F_INVALID, EDOM, x, 0.0F,1);
 71            else
 72              /* Return positive infinity with no signal */
 73              return x;
 74          }
 75      }
 76    else if ((ux & SIGNBIT_SP32) || (ux < 0x3f800000))
 77      {
 78        /* x is less than 1.0. Return a NaN. */
 79        return __amd_handle_errorf(_FUNCNAME,__amd_acosh, INDEFBITPATT_SP32, _DOMAIN,
 80                             AMD_F_INVALID, EDOM, x, 0.0F,1);
 81      }
 82  
 83    dx = x;
 84  
 85    if (ux > recrteps)
 86      {
 87        /* Arguments greater than 1/sqrt(epsilon) in magnitude are
 88           approximated by acoshf(x) = ln(2) + ln(x) */
 89        r = FN_PROTOTYPE(log)(dx) + log2;
 90      }
 91    else if (ux > 0x40000000)
 92      {
 93        /* 2.0 <= x <= 1/sqrt(epsilon) */
 94        /* acoshf for these arguments is approximated by
 95           acoshf(x) = ln(x + sqrt(x*x-1)) */
 96        rarg = dx*dx-1.0;
 97        /* Use assembly instruction to compute r = sqrt(rarg); */
 98        ASMSQRT(rarg,r);
 99        rarg = r + dx;
100        r = FN_PROTOTYPE(log)(rarg);
101      }
102    else
103      {
104        /* sqrt(epsilon) <= x <= 2.0 */
105        t = dx - 1.0;
106        rarg = 2.0*t + t*t;
107        ASMSQRT(rarg,r);  /* r = sqrt(rarg) */
108        rarg = t + r;
109        r = FN_PROTOTYPE(log1p)(rarg);
110      }
111    return (float)(r);
112  }
113  
114