
/*
*  Copyright (C) 2008-2009 Advanced Micro Devices, Inc. All Rights Reserved.
*
*  This file is part of libacml_mv.
*
*  libacml_mv is free software; you can redistribute it and/or
*  modify it under the terms of the GNU Lesser General Public
*  License as published by the Free Software Foundation; either
*  version 2.1 of the License, or (at your option) any later version.
*
*  libacml_mv is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
*  Lesser General Public License for more details.
*
*  You should have received a copy of the GNU Lesser General Public
*  License along with libacml_mv.  If not, see
*  <http://www.gnu.org/licenses/>.
*
*/


#include "../inc/libm_amd.h"
#include "../inc/libm_util_amd.h"

#define USE_HANDLE_ERROR
#define USE_LOG_KERNEL_AMD
#define USE_VAL_WITH_FLAGS
#include "../inc/libm_inlines_amd.h"
#undef USE_HANDLE_ERROR
#undef USE_LOG_KERNEL_AMD
#undef VAL_WITH_FLAGS

#undef _FUNCNAME
#define _FUNCNAME "asinh"
double FN_PROTOTYPE(asinh)(double x)
{

  unsigned long long ux, ax, xneg;
  double absx, r, rarg, t, r1, r2, poly, s, v1, v2;
  int xexp;

  static const unsigned long long
    rteps = 0x3e46a09e667f3bcd,    /* sqrt(eps) = 1.05367121277235086670e-08 */
    recrteps = 0x4196a09e667f3bcd; /* 1/rteps = 9.49062656242515593767e+07 */

  /* log2_lead and log2_tail sum to an extra-precise version
     of log(2) */
  static const double
    log2_lead = 6.93147122859954833984e-01,  /* 0x3fe62e42e0000000 */
    log2_tail = 5.76999904754328540596e-08;  /* 0x3e6efa39ef35793c */


  GET_BITS_DP64(x, ux);
  ax = ux & ~SIGNBIT_DP64;
  xneg = ux & SIGNBIT_DP64;
  PUT_BITS_DP64(ax, absx);

  if ((ux & EXPBITS_DP64) == EXPBITS_DP64)
    {
      /* x is either NaN or infinity */
      if (ux & MANTBITS_DP64)
        {
          /* x is NaN */
#ifdef WINDOWS
          return handle_error(_FUNCNAME, ux|0x0008000000000000, _DOMAIN,
                              AMD_F_INVALID, EDOM, x, 0.0);
#else
          return x + x; /* Raise invalid if it is a signalling NaN */
#endif
        }
      else
        {
          /* x is infinity. Return the same infinity. */
#ifdef WINDOWS
          if (ux & SIGNBIT_DP64)
            return handle_error(_FUNCNAME, NINFBITPATT_DP64, _DOMAIN,
                                AMD_F_INVALID, EDOM, x, 0.0);
          else
            return handle_error(_FUNCNAME, PINFBITPATT_DP64, _DOMAIN,
                                AMD_F_INVALID, EDOM, x, 0.0);
#else
          return x;
#endif
        }
    }
  else if (ax < rteps) /* abs(x) < sqrt(epsilon) */
    {
      if (ax == 0x0000000000000000)
        {
          /* x is +/-zero. Return the same zero. */
          return x;
        }
      else
        {
          /* Tiny arguments approximated by asinh(x) = x
             - avoid slow operations on denormalized numbers */
          return val_with_flags(x,AMD_F_INEXACT);
        }
    }


  if (ax <= 0x3ff0000000000000) /* abs(x) <= 1.0 */
    {
      /* Arguments less than 1.0 in magnitude are
         approximated by [4,4] or [5,4] minimax polynomials
         fitted to asinh series 4.6.31 (x < 1) from Abramowitz and Stegun
      */
      t = x*x;
      if (ax < 0x3fd0000000000000)
        {
          /* [4,4] for 0 < abs(x) < 0.25 */
          poly =
            (-0.12845379283524906084997e0 +
             (-0.21060688498409799700819e0 +
              (-0.10188951822578188309186e0 +
               (-0.13891765817243625541799e-1 -
                0.10324604871728082428024e-3 * t) * t) * t) * t) /
            (0.77072275701149440164511e0 +
             (0.16104665505597338100747e1 +
              (0.11296034614816689554875e1 +
               (0.30079351943799465092429e0 +
                0.235224464765951442265117e-1 * t) * t) * t) * t);
        }
      else if (ax < 0x3fe0000000000000)
        {
          /* [4,4] for 0.25 <= abs(x) < 0.5 */
          poly =
            (-0.12186605129448852495563e0 +
             (-0.19777978436593069928318e0 +
              (-0.94379072395062374824320e-1 +
               (-0.12620141363821680162036e-1 -
                0.903396794842691998748349e-4 * t) * t) * t) * t) /
            (0.73119630776696495279434e0 +
             (0.15157170446881616648338e1 +
              (0.10524909506981282725413e1 +
               (0.27663713103600182193817e0 +
                0.21263492900663656707646e-1 * t) * t) * t) * t);
        }
      else if (ax < 0x3fe8000000000000)
        {
          /* [4,4] for 0.5 <= abs(x) < 0.75 */
          poly =
            (-0.81210026327726247622500e-1 +
             (-0.12327355080668808750232e0 +
              (-0.53704925162784720405664e-1 +
               (-0.63106739048128554465450e-2 -
                0.35326896180771371053534e-4 * t) * t) * t) * t) /
            (0.48726015805581794231182e0 +
             (0.95890837357081041150936e0 +
              (0.62322223426940387752480e0 +
               (0.15028684818508081155141e0 +
                0.10302171620320141529445e-1 * t) * t) * t) * t);
        }
      else
        {
          /* [5,4] for 0.75 <= abs(x) <= 1.0 */
          poly =
            (-0.4638179204422665073e-1 +
             (-0.7162729496035415183e-1 +
              (-0.3247795155696775148e-1 +
               (-0.4225785421291932164e-2 +
                (-0.3808984717603160127e-4 +
                 0.8023464184964125826e-6 * t) * t) * t) * t) * t) /
            (0.2782907534642231184e0 +
             (0.5549945896829343308e0 +
              (0.3700732511330698879e0 +
               (0.9395783438240780722e-1 +
                0.7200057974217143034e-2 * t) * t) * t) * t);
        }
      return x + x*t*poly;
    }
  else if (ax < 0x4040000000000000)
    {
      /* 1.0 <= abs(x) <= 32.0 */
      /* Arguments in this region are approximated by various
         minimax polynomials fitted to asinh series 4.6.31
         in Abramowitz and Stegun.
      */
      t = x*x;
      if (ax >= 0x4020000000000000)
        {
          /* [3,3] for 8.0 <= abs(x) <= 32.0 */
          poly =
            (-0.538003743384069117e-10 +
             (-0.273698654196756169e-9 +
              (-0.268129826956403568e-9 -
               0.804163374628432850e-29 * t) * t) * t) /
            (0.238083376363471960e-9 +
             (0.203579344621125934e-8 +
              (0.450836980450693209e-8 +
               0.286005148753497156e-8 * t) * t) * t);
        }
      else if (ax >= 0x4010000000000000)
        {
          /* [4,3] for 4.0 <= abs(x) <= 8.0 */
          poly =
            (-0.178284193496441400e-6 +
             (-0.928734186616614974e-6 +
              (-0.923318925566302615e-6 +
               (-0.776417026702577552e-19 +
                0.290845644810826014e-21 * t) * t) * t) * t) /
            (0.786694697277890964e-6 +
             (0.685435665630965488e-5 +
              (0.153780175436788329e-4 +
               0.984873520613417917e-5 * t) * t) * t);

        }
      else if (ax >= 0x4000000000000000)
        {
          /* [5,4] for 2.0 <= abs(x) <= 4.0 */
          poly =
            (-0.209689451648100728e-6 +
             (-0.219252358028695992e-5 +
              (-0.551641756327550939e-5 +
               (-0.382300259826830258e-5 +
                (-0.421182121910667329e-17 +
                 0.492236019998237684e-19 * t) * t) * t) * t) * t) /
            (0.889178444424237735e-6 +
             (0.131152171690011152e-4 +
              (0.537955850185616847e-4 +
               (0.814966175170941864e-4 +
                0.407786943832260752e-4 * t) * t) * t) * t);
        }
      else if (ax >= 0x3ff8000000000000)
        {
          /* [5,4] for 1.5 <= abs(x) <= 2.0 */
          poly =
            (-0.195436610112717345e-4 +
             (-0.233315515113382977e-3 +
              (-0.645380957611087587e-3 +
               (-0.478948863920281252e-3 +
                (-0.805234112224091742e-12 +
                 0.246428598194879283e-13 * t) * t) * t) * t) * t) /
            (0.822166621698664729e-4 +
             (0.135346265620413852e-2 +
              (0.602739242861830658e-2 +
               (0.972227795510722956e-2 +
                0.510878800983771167e-2 * t) * t) * t) * t);
        }
      else
        {
          /* [5,5] for 1.0 <= abs(x) <= 1.5 */
          poly =
            (-0.121224194072430701e-4 +
             (-0.273145455834305218e-3 +
              (-0.152866982560895737e-2 +
               (-0.292231744584913045e-2 +
                (-0.174670900236060220e-2 -
                 0.891754209521081538e-12 * t) * t) * t) * t) * t) /
            (0.499426632161317606e-4 +
             (0.139591210395547054e-2 +
              (0.107665231109108629e-1 +
               (0.325809818749873406e-1 +
                (0.415222526655158363e-1 +
                 0.186315628774716763e-1 * t) * t) * t) * t) * t);
        }
      log_kernel_amd64(absx, ax, &xexp, &r1, &r2);
      r1 = ((xexp+1) * log2_lead + r1);
      r2 = ((xexp+1) * log2_tail + r2);
      /* Now (r1,r2) sum to log(2x). Add the term
         1/(2.2.x^2) = 0.25/t, and add poly/t, carefully
         to maintain precision. (Note that we add poly/t
         rather than poly because of the *x factor used
         when generating the minimax polynomial) */
      v2 = (poly+0.25)/t;
      r = v2 + r1;
      s = ((r1 - r) + v2) + r2;
      v1 = r + s;
      v2 = (r - v1) + s;
      r = v1 + v2;
      if (xneg)
        return -r;
      else
        return r;
    }
  else
    {
      /* abs(x) > 32.0 */
      if (ax > recrteps)
        {
          /* Arguments greater than 1/sqrt(epsilon) in magnitude are
             approximated by asinh(x) = ln(2) + ln(abs(x)), with sign of x */
          /* log_kernel_amd(x) returns xexp, r1, r2 such that
             log(x) = xexp*log(2) + r1 + r2 */
          log_kernel_amd64(absx, ax, &xexp, &r1, &r2);
          /* Add (xexp+1) * log(2) to z1,z2 to get the result asinh(x).
             The computed r1 is not subject to rounding error because
             (xexp+1) has at most 10 significant bits, log(2) has 24 significant
             bits, and r1 has up to 24 bits; and the exponents of r1
             and r2 differ by at most 6. */
          r1 = ((xexp+1) * log2_lead + r1);
          r2 = ((xexp+1) * log2_tail + r2);
          if (xneg)
            return -(r1 + r2);
          else
            return r1 + r2;
        }
      else
        {
          rarg = absx*absx+1.0;
          /* Arguments such that 32.0 <= abs(x) <= 1/sqrt(epsilon) are
             approximated by
               asinh(x) = ln(abs(x) + sqrt(x*x+1))
             with the sign of x (see Abramowitz and Stegun 4.6.20) */
          /* Use assembly instruction to compute r = sqrt(rarg); */
          ASMSQRT(rarg,r);
          r += absx;
          GET_BITS_DP64(r, ax);
          log_kernel_amd64(r, ax, &xexp, &r1, &r2);
          r1 = (xexp * log2_lead + r1);
          r2 = (xexp * log2_tail + r2);
          if (xneg)
            return -(r1 + r2);
          else
            return r1 + r2;
        }
    }
}

weak_alias (__asinh, asinh)