src/atanhf.c

/*
*  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"

#include <stdio.h>

#define USE_NANF_WITH_FLAGS
#define USE_VALF_WITH_FLAGS
#define USE_INFINITYF_WITH_FLAGS
#define USE_HANDLE_ERRORF
#include "../inc/libm_inlines_amd.h"
#undef USE_NANF_WITH_FLAGS
#undef USE_VALF_WITH_FLAGS
#undef USE_INFINITYF_WITH_FLAGS
#undef USE_HANDLE_ERRORF

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

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

#undef _FUNCNAME
#define _FUNCNAME "atanhf"
float FN_PROTOTYPE(atanhf)(float x)
{

  double dx;
  unsigned int ux, ax;
  double r, t, poly;

  GET_BITS_SP32(x, ux);
  ax = ux & ~SIGNBIT_SP32;

  if ((ux & EXPBITS_SP32) == EXPBITS_SP32)
    {
      /* x is either NaN or infinity */
      if (ux & MANTBITS_SP32)
        {
          /* x is NaN */
#ifdef WINDOWS
          return handle_errorf(_FUNCNAME, ux|0x00400000, _DOMAIN,
                              0, EDOM, x, 0.0F);
#else
          return x + x; /* Raise invalid if it is a signalling NaN */
#endif
        }
      else
        {
          /* x is infinity; return a NaN */
#ifdef WINDOWS
          return handle_errorf(_FUNCNAME, INDEFBITPATT_SP32, _DOMAIN,
                               AMD_F_INVALID, EDOM, x, 0.0F);
#else
          return retval_errno_edom(x,nanf_with_flags(AMD_F_INVALID));
#endif
        }
    }
  else if (ax >= 0x3f800000)
    {
      if (ax > 0x3f800000)
        {
          /* abs(x) > 1.0; return NaN */
#ifdef WINDOWS
          return handle_errorf(_FUNCNAME, INDEFBITPATT_SP32, _DOMAIN,
                               AMD_F_INVALID, EDOM, x, 0.0F);
#else
          return retval_errno_edom(x,nanf_with_flags(AMD_F_INVALID));
#endif
        }
      else if (ux == 0x3f800000)
        {
          /* x = +1.0; return infinity with the same sign as x
             and set the divbyzero status flag */
#ifdef WINDOWS
          return handle_errorf(_FUNCNAME, PINFBITPATT_SP32, _DOMAIN,
                               AMD_F_INVALID, EDOM, x, 0.0F);
#else
          return retval_errno_edom(x,infinityf_with_flags(AMD_F_DIVBYZERO));
#endif
        }
      else
        {
          /* x = -1.0; return infinity with the same sign as x */
#ifdef WINDOWS
          return handle_errorf(_FUNCNAME, NINFBITPATT_SP32, _DOMAIN,
                               AMD_F_INVALID, EDOM, x, 0.0F);
#else
          return retval_errno_edom(x,-infinityf_with_flags(AMD_F_DIVBYZERO));
#endif
        }
    }

  if (ax < 0x39000000)
    {
      if (ax == 0x00000000)
        {
          /* x is +/-zero. Return the same zero. */
          return x;
        }
      else
        {
          /* Arguments smaller than 2^(-13) in magnitude are
             approximated by atanhf(x) = x, raising inexact flag. */
          return valf_with_flags(x, AMD_F_INEXACT);
        }
    }
  else
    {
      dx = x;
      if (ax < 0x3f000000)
        {
          /* Arguments up to 0.5 in magnitude are
             approximated by a [2,2] minimax polynomial */
          t = dx*dx;
          poly =
            (0.39453629046e0 +
           (-0.28120347286e0 +
             0.92834212715e-2 * t) * t) /
            (0.11836088638e1 + 
           (-0.15537744551e1 +
             0.45281890445e0 * t) * t);
          return (float)(dx + dx*t*poly);
        }
      else
        {
          /* abs(x) >= 0.5 */
          /* Note that
               atanhf(x) = 0.5 * ln((1+x)/(1-x))
             (see Abramowitz and Stegun 4.6.22).
             For greater accuracy we use the variant formula
             atanhf(x) = log(1 + 2x/(1-x)) = log1p(2x/(1-x)).
          */
          if (ux & SIGNBIT_SP32)
            {
              /* Argument x is negative */
              r = (-2.0 * dx) / (1.0 + dx);
              r = 0.5 * FN_PROTOTYPE(log1p)(r);
              return (float)-r;
            }
          else
            {
              r = (2.0 * dx) / (1.0 - dx);
              r = 0.5 * FN_PROTOTYPE(log1p)(r);
              return (float)r;
            }
        }
    }
}

weak_alias (__atanhf, atanhf)