blob: 38692b4454151cc1aef454f83795d9f89e1cb9a9 [file] [log] [blame]
/*
* 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)