
/*
*  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_NAN_WITH_FLAGS
#define USE_HANDLE_ERROR
#define USE_LOG_KERNEL_AMD
#include "../inc/libm_inlines_amd.h"
#undef USE_NAN_WITH_FLAGS
#undef USE_HANDLE_ERROR
#undef USE_LOG_KERNEL_AMD

#include "../inc/libm_errno_amd.h"

#ifndef WINDOWS
/* Deal with errno for out-of-range argument */
static inline double retval_errno_edom(double x)
{
  struct exception exc;
  exc.arg1 = x;
  exc.arg2 = x;
  exc.type = DOMAIN;
  exc.name = (char *)"acosh";
  if (_LIB_VERSION == _SVID_)
    exc.retval = -HUGE;
  else
    exc.retval = nan_with_flags(AMD_F_INVALID);
  if (_LIB_VERSION == _POSIX_)
    __set_errno(EDOM);
  else if (!matherr(&exc))
    {
      if(_LIB_VERSION == _SVID_)
        (void)fputs("acosh: DOMAIN error\n", stderr);
    __set_errno(EDOM);
    }
  return exc.retval;
}
#endif

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

  unsigned long long ux;
  double r, rarg, r1, r2;
  int xexp;

  static const unsigned long long
    recrteps = 0x4196a09e667f3bcd; /* 1/sqrt(eps) = 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);

  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 */
          if (ux & SIGNBIT_DP64)
            /* x is negative infinity. Return a NaN. */
#ifdef WINDOWS
            return handle_error(_FUNCNAME, INDEFBITPATT_DP64, _DOMAIN,
                                AMD_F_INVALID, EDOM, x, 0.0);
#else
            return retval_errno_edom(x);
#endif
          else
            /* Return positive infinity with no signal */
            return x;
        }
    }
  else if ((ux & SIGNBIT_DP64) || (ux <= 0x3ff0000000000000))
    {
      /* x <= 1.0 */
      if (ux == 0x3ff0000000000000)
        {
          /* x = 1.0; return zero. */
          return 0.0;
        }
      else
        {
          /* x is less than 1.0. Return a NaN. */
#ifdef WINDOWS
          return handle_error(_FUNCNAME, INDEFBITPATT_DP64, _DOMAIN,
                              AMD_F_INVALID, EDOM, x, 0.0);
#else
          return retval_errno_edom(x);
#endif
        }
    }


  if (ux > recrteps)
    {
      /* Arguments greater than 1/sqrt(epsilon) in magnitude are
         approximated by acosh(x) = ln(2) + ln(x) */
      /* log_kernel_amd(x) returns xexp, r1, r2 such that
         log(x) = xexp*log(2) + r1 + r2 */
      log_kernel_amd64(x, ux, &xexp, &r1, &r2);
      /* Add (xexp+1) * log(2) to z1,z2 to get the result acosh(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);
      return r1 + r2;
    }
  else if (ux >= 0x4060000000000000)
    {
      /* 128.0 <= x <= 1/sqrt(epsilon) */
      /* acosh for these arguments is approximated by
         acosh(x) = ln(x + sqrt(x*x-1)) */
      rarg = x*x-1.0;
      /* Use assembly instruction to compute r = sqrt(rarg); */
      ASMSQRT(rarg,r);
      r += x;
      GET_BITS_DP64(r, ux);
      log_kernel_amd64(r, ux, &xexp, &r1, &r2);
      r1 = (xexp * log2_lead + r1);
      r2 = (xexp * log2_tail + r2);
      return r1 + r2;
    }
  else
    {
      /* 1.0 < x <= 128.0 */
      double u1, u2, v1, v2, w1, w2, hx, tx, t, r, s, p1, p2, a1, a2, c1, c2,
        poly;
      if (ux >= 0x3ff8000000000000)
        {
          /* 1.5 <= x <= 128.0 */
          /* We use minimax polynomials,
             based on Abramowitz and Stegun 4.6.32 series
             expansion for acosh(x), with the log(2x) and 1/(2.2.x^2)
             terms removed. We compensate for these two terms later.
          */
          t = x*x;
          if (ux >= 0x4040000000000000)
            {
              /* [3,2] for 32.0 <= x <= 128.0 */
              poly =
                (0.45995704464157438175e-9 +
                 (-0.89080839823528631030e-9 +
                  (-0.10370522395596168095e-27 +
                   0.35255386405811106347e-32 * t) * t) * t) /
                (0.21941191335882074014e-8 +
                 (-0.10185073058358334569e-7 +
                  0.95019562478430648685e-8 * t) * t);
            }
          else if (ux >= 0x4020000000000000)
            {
              /* [3,3] for 8.0 <= x <= 32.0 */
              poly =
                (-0.54903656589072526589e-10 +
                 (0.27646792387218569776e-9 +
                  (-0.26912957240626571979e-9 -
                   0.86712268396736384286e-29 * t) * t) * t) /
                (-0.24327683788655520643e-9 +
                 (0.20633757212593175571e-8 +
                  (-0.45438330985257552677e-8 +
                   0.28707154390001678580e-8 * t) * t) * t);
            }
          else if (ux >= 0x4010000000000000)
            {
              /* [4,3] for 4.0 <= x <= 8.0 */
              poly =
                (-0.20827370596738166108e-6 +
                 (0.10232136919220422622e-5 +
                  (-0.98094503424623656701e-6 +
                   (-0.11615338819596146799e-18 +
                    0.44511847799282297160e-21 * t) * t) * t) * t) /
                (-0.92579451630913718588e-6 +
                 (0.76997374707496606639e-5 +
                  (-0.16727286999128481170e-4 +
                   0.10463413698762590251e-4 * t) * t) * t);
            }
          else if (ux >= 0x4000000000000000)
            {
              /* [5,5] for 2.0 <= x <= 4.0 */
              poly =
                (-0.122195030526902362060e-7 +
                 (0.157894522814328933143e-6 +
                  (-0.579951798420930466109e-6 +
                   (0.803568881125803647331e-6 +
                    (-0.373906657221148667374e-6 -
                     0.317856399083678204443e-21 * t) * t) * t) * t) * t) /
                (-0.516260096352477148831e-7 +
                 (0.894662592315345689981e-6 +
                  (-0.475662774453078218581e-5 +
                   (0.107249291567405130310e-4 +
                    (-0.107871445525891289759e-4 +
                     0.398833767702587224253e-5 * t) * t) * t) * t) * t);
            }
          else if (ux >= 0x3ffc000000000000)
            {
              /* [5,4] for 1.75 <= x <= 2.0 */
              poly =
                (0.1437926821253825186e-3 +
                 (-0.1034078230246627213e-2 +
                  (0.2015310005461823437e-2 +
                   (-0.1159685218876828075e-2 +
                    (-0.9267353551307245327e-11 +
                     0.2880267770324388034e-12 * t) * t) * t) * t) * t) /
                (0.6305521447028109891e-3 +
                 (-0.6816525887775002944e-2 +
                  (0.2228081831550003651e-1 +
                   (-0.2836886105406603318e-1 +
                    0.1236997707206036752e-1 * t) * t) * t) * t);
            }
          else
            {
              /* [5,4] for 1.5 <= x <= 1.75 */
              poly =
                ( 0.7471936607751750826e-3 +
                  (-0.4849405284371905506e-2 +
                   (0.8823068059778393019e-2 +
                    (-0.4825395461288629075e-2 +
                     (-0.1001984320956564344e-8 +
                      0.4299919281586749374e-10 * t) * t) * t) * t) * t) /
                (0.3322359141239411478e-2 +
                 (-0.3293525930397077675e-1 +
                  (0.1011351440424239210e0 +
                   (-0.1227083591622587079e0 +
                    0.5147099404383426080e-1 * t) * t) * t) * t);
            }
          GET_BITS_DP64(x, ux);
          log_kernel_amd64(x, ux, &xexp, &r1, &r2);
          r1 = ((xexp+1) * log2_lead + r1);
          r2 = ((xexp+1) * log2_tail + r2);
          /* Now (r1,r2) sum to log(2x). Subtract 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;
          return v1 + ((r - v1) + s);
        }

      /* Here 1.0 <= x <= 1.5. It is hard to maintain accuracy here so
         we have to go to great lengths to do so. */

      /* We compute the value
           t = x - 1.0 + sqrt(2.0*(x - 1.0) + (x - 1.0)*(x - 1.0))
         using simulated quad precision. */
      t = x - 1.0;
      u1 = t * 2.0;

      /* dekker_mul12(t,t,&v1,&v2); */
      GET_BITS_DP64(t, ux);
      ux &= 0xfffffffff8000000;
      PUT_BITS_DP64(ux, hx);
      tx = t - hx;
      v1 = t * t;
      v2 = (((hx * hx - v1) + hx * tx) + tx * hx) + tx * tx;

      /* dekker_add2(u1,0.0,v1,v2,&w1,&w2); */
      r = u1 + v1;
      s = (((u1 - r) + v1) + v2);
      w1 = r + s;
      w2 = (r - w1) + s;

      /* dekker_sqrt2(w1,w2,&u1,&u2); */
      ASMSQRT(w1,p1);
      GET_BITS_DP64(p1, ux);
      ux &= 0xfffffffff8000000;
      PUT_BITS_DP64(ux, c1);
      c2 = p1 - c1;
      a1 = p1 * p1;
      a2 = (((c1 * c1 - a1) + c1 * c2) + c2 * c1) + c2 * c2;
      p2 = (((w1 - a1) - a2) + w2) * 0.5 / p1;
      u1 = p1 + p2;
      u2 = (p1 - u1) + p2;

      /* dekker_add2(u1,u2,t,0.0,&v1,&v2); */
      r = u1 + t;
      s = (((u1 - r) + t)) + u2;
      r1 = r + s;
      r2 = (r - r1) + s;
      t = r1 + r2;

      /* Check for x close to 1.0. */
      if (x < 1.13)
        {
          /* Here 1.0 <= x < 1.13 implies r <= 0.656. In this region
             we need to take extra care to maintain precision.
             We have t = r1 + r2 = (x - 1.0 + sqrt(x*x-1.0))
             to more than basic precision. We use the Taylor series
             for log(1+x), with terms after the O(x*x) term
             approximated by a [6,6] minimax polynomial. */
          double b1, b2, c1, c2, e1, e2, q1, q2, c, cc, hr1, tr1, hpoly, tpoly, hq1, tq1, hr2, tr2;
          poly =
            (0.30893760556597282162e-21 +
             (0.10513858797132174471e0 +
              (0.27834538302122012381e0 +
               (0.27223638654807468186e0 +
                (0.12038958198848174570e0 +
                 (0.23357202004546870613e-1 +
                  (0.15208417992520237648e-2 +
                   0.72741030690878441996e-7 * t) * t) * t) * t) * t) * t) * t) /
            (0.31541576391396523486e0 +
             (0.10715979719991342022e1 +
              (0.14311581802952004012e1 +
               (0.94928647994421895988e0 +
                (0.32396235926176348977e0 +
                 (0.52566134756985833588e-1 +
                  0.30477895574211444963e-2 * t) * t) * t) * t) * t) * t);

          /* Now we can compute the result r = acosh(x) = log1p(t)
             using the formula t - 0.5*t*t + poly*t*t. Since t is
             represented as r1+r2, the formula becomes
             r = r1+r2 - 0.5*(r1+r2)*(r1+r2) + poly*(r1+r2)*(r1+r2).
             Expanding out, we get
               r = r1 + r2 - (0.5 + poly)*(r1*r1 + 2*r1*r2 + r2*r2)
             and ignoring negligible quantities we get
               r = r1 + r2 - 0.5*r1*r1 + r1*r2 + poly*t*t
          */
          if (x < 1.06)
            {
              double b, c, e;
              b = r1*r2;
              c = 0.5*r1*r1;
              e = poly*t*t;
              /* N.B. the order of additions and subtractions is important */
              r = (((r2 - b) + e) - c) + r1;
              return r;
            }
          else
            {
              /* For 1.06 <= x <= 1.13 we must evaluate in extended precision
                 to reach about 1 ulp accuracy (in this range the simple code
                 above only manages about 1.5 ulp accuracy) */

              /* Split poly, r1 and r2 into head and tail sections */
              GET_BITS_DP64(poly, ux);
              ux &= 0xfffffffff8000000;
              PUT_BITS_DP64(ux,hpoly);
              tpoly = poly - hpoly;
              GET_BITS_DP64(r1,ux);
              ux &= 0xfffffffff8000000;
              PUT_BITS_DP64(ux,hr1);
              tr1 = r1 - hr1;
              GET_BITS_DP64(r2, ux);
              ux &= 0xfffffffff8000000;
              PUT_BITS_DP64(ux,hr2);
              tr2 = r2 - hr2;

              /* e = poly*t*t */
              c = poly * r1;
              cc = (((hpoly * hr1 - c) + hpoly * tr1) + tpoly * hr1) + tpoly * tr1;
              cc = poly * r2 + cc;
              q1 = c + cc;
              q2 = (c - q1) + cc;
              GET_BITS_DP64(q1, ux);
              ux &= 0xfffffffff8000000;
              PUT_BITS_DP64(ux,hq1);
              tq1 = q1 - hq1;
              c = q1 * r1;
              cc = (((hq1 * hr1 - c) + hq1 * tr1) + tq1 * hr1) + tq1 * tr1;
              cc = q1 * r2 + q2 * r1 + cc;
              e1 = c + cc;
              e2 = (c - e1) + cc;

              /* b = r1*r2 */
              b1 = r1 * r2;
              b2 = (((hr1 * hr2 - b1) + hr1 * tr2) + tr1 * hr2) + tr1 * tr2;

              /* c = 0.5*r1*r1 */
              c1 = (0.5*r1) * r1;
              c2 = (((0.5*hr1 * hr1 - c1) + 0.5*hr1 * tr1) + 0.5*tr1 * hr1) + 0.5*tr1 * tr1;

              /* v = a + d - b */
              r = r1 - b1;
              s = (((r1 - r) - b1) - b2) + r2;
              v1 = r + s;
              v2 = (r - v1) + s;

              /* w = (a + d - b) - c */
              r = v1 - c1;
              s = (((v1 - r) - c1) - c2) + v2;
              w1 = r + s;
              w2 = (r - w1) + s;

              /* u = ((a + d - b) - c) + e */
              r = w1 + e1;
              s = (((w1 - r) + e1) + e2) + w2;
              u1 = r + s;
              u2 = (r - u1) + s;

              /* The result r = acosh(x) */
              r = u1 + u2;

              return r;
            }
        }
      else
        {
          /* For arguments 1.13 <= x <= 1.5 the log1p function
             is good enough */
            return FN_PROTOTYPE(log1p)(t);
        }
    }
}

weak_alias (__acosh, acosh)