lib/libm/src/e_atanhl.c

/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */

/* @(#)e_atanh.c 1.3 95/01/18 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 *
 * $FreeBSD: head/lib/msun/src/e_atanhl.c 251599 2013-06-10 06:04:58Z das $
 */

/*
 * See e_atanh.c for complete comments.
 *
 * Converted to long double by David Schultz <das@FreeBSD.ORG> and
 * Bruce D. Evans.
 */

#include <float.h>
#ifdef __i386__
#include <ieeefp.h>
#endif

#include "fpmath.h"
#include "math.h"
#include "math_private.h"

/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
#if LDBL_MANT_DIG == 64
#define	EXP_TINY	-34
#elif LDBL_MANT_DIG == 113
#define	EXP_TINY	-58
#else
#error "Unsupported long double format"
#endif

#if LDBL_MAX_EXP != 0x4000
/* We also require the usual expsign encoding. */
#error "Unsupported long double format"
#endif

#define	BIAS	(LDBL_MAX_EXP - 1)

static const double one = 1.0, huge = 1e300;
static const double zero = 0.0;

long double
atanhl(long double x)
{
	long double t;
	uint16_t hx, ix;

	ENTERI();
	GET_LDBL_EXPSIGN(hx, x);
	ix = hx & 0x7fff;
	if (ix >= 0x3fff)		/* |x| >= 1, or NaN or misnormal */
	    RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
	if (ix < BIAS + EXP_TINY && (huge + x) > zero)
	    RETURNI(x);			/* x is tiny */
	SET_LDBL_EXPSIGN(x, ix);
	if (ix < 0x3ffe) {		/* |x| < 0.5, or misnormal */
	    t = x+x;
	    t = 0.5*log1pl(t+t*x/(one-x));
	} else
	    t = 0.5*log1pl((x+x)/(one-x));
	RETURNI((hx & 0x8000) == 0 ? t : -t);
}
Commit	Line	Data
	1	/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
	2
	3	/* @(#)e_atanh.c 1.3 95/01/18 */
	4	/*
	5	* ====================================================
	6	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
	7	*
	8	* Developed at SunSoft, a Sun Microsystems, Inc. business.
	9	* Permission to use, copy, modify, and distribute this
	10	* software is freely granted, provided that this notice
	11	* is preserved.
	12	* ====================================================
	13	*
	14	* $FreeBSD: head/lib/msun/src/e_atanhl.c 251599 2013-06-10 06:04:58Z das $
	15	*/
	16
	17	/*
	18	* See e_atanh.c for complete comments.
	19	*
	20	* Converted to long double by David Schultz <das@FreeBSD.ORG> and
	21	* Bruce D. Evans.
	22	*/
	23
	24	#include <float.h>
	25	#ifdef __i386__
	26	#include <ieeefp.h>
	27	#endif
	28
	29	#include "fpmath.h"
	30	#include "math.h"
	31	#include "math_private.h"
	32
	33	/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
	34	#if LDBL_MANT_DIG == 64
	35	#define EXP_TINY -34
	36	#elif LDBL_MANT_DIG == 113
	37	#define EXP_TINY -58
	38	#else
	39	#error "Unsupported long double format"
	40	#endif
	41
	42	#if LDBL_MAX_EXP != 0x4000
	43	/* We also require the usual expsign encoding. */
	44	#error "Unsupported long double format"
	45	#endif
	46
	47	#define BIAS (LDBL_MAX_EXP - 1)
	48
	49	static const double one = 1.0, huge = 1e300;
	50	static const double zero = 0.0;
	51
	52	long double
	53	atanhl(long double x)
	54	{
	55	long double t;
	56	uint16_t hx, ix;
	57
	58	ENTERI();
	59	GET_LDBL_EXPSIGN(hx, x);
	60	ix = hx & 0x7fff;
	61	if (ix >= 0x3fff) /* \|x\| >= 1, or NaN or misnormal */
	62	RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
	63	if (ix < BIAS + EXP_TINY && (huge + x) > zero)
	64	RETURNI(x); /* x is tiny */
	65	SET_LDBL_EXPSIGN(x, ix);
	66	if (ix < 0x3ffe) { /* \|x\| < 0.5, or misnormal */
	67	t = x+x;
	68	t = 0.5log1pl(t+tx/(one-x));
	69	} else
	70	t = 0.5*log1pl((x+x)/(one-x));
	71	RETURNI((hx & 0x8000) == 0 ? t : -t);
	72	}