1251599Sdas/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
2141296Sdas
3141296Sdas/* @(#)e_atanh.c 1.3 95/01/18 */
42116Sjkh/*
52116Sjkh * ====================================================
62116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
72116Sjkh *
8141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business.
92116Sjkh * Permission to use, copy, modify, and distribute this
10141296Sdas * software is freely granted, provided that this notice
112116Sjkh * is preserved.
122116Sjkh * ====================================================
13141296Sdas *
142116Sjkh */
152116Sjkh
16176451Sdas#include <sys/cdefs.h>
17176451Sdas__FBSDID("$FreeBSD$");
182116Sjkh
19251599Sdas/*
20251599Sdas * See e_atanh.c for complete comments.
212116Sjkh *
22251599Sdas * Converted to long double by David Schultz <das@FreeBSD.ORG> and
23251599Sdas * Bruce D. Evans.
242116Sjkh */
252116Sjkh
26251599Sdas#include <float.h>
27251599Sdas#ifdef __i386__
28251599Sdas#include <ieeefp.h>
29251599Sdas#endif
30251599Sdas
31251599Sdas#include "fpmath.h"
322116Sjkh#include "math.h"
332116Sjkh#include "math_private.h"
342116Sjkh
35251599Sdas/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
36251599Sdas#if LDBL_MANT_DIG == 64
37251599Sdas#define	EXP_TINY	-34
38251599Sdas#elif LDBL_MANT_DIG == 113
39251599Sdas#define	EXP_TINY	-58
40251599Sdas#else
41251599Sdas#error "Unsupported long double format"
42251599Sdas#endif
43251599Sdas
44251599Sdas#if LDBL_MAX_EXP != 0x4000
45251599Sdas/* We also require the usual expsign encoding. */
46251599Sdas#error "Unsupported long double format"
47251599Sdas#endif
48251599Sdas
49251599Sdas#define	BIAS	(LDBL_MAX_EXP - 1)
50251599Sdas
512116Sjkhstatic const double one = 1.0, huge = 1e300;
522116Sjkhstatic const double zero = 0.0;
532116Sjkh
54251599Sdaslong double
55251599Sdasatanhl(long double x)
562116Sjkh{
57251599Sdas	long double t;
58251599Sdas	uint16_t hx, ix;
59251599Sdas
60251599Sdas	ENTERI();
61251599Sdas	GET_LDBL_EXPSIGN(hx, x);
62251599Sdas	ix = hx & 0x7fff;
63251599Sdas	if (ix >= 0x3fff)		/* |x| >= 1, or NaN or misnormal */
64251599Sdas	    RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
65251599Sdas	if (ix < BIAS + EXP_TINY && (huge + x) > zero)
66251599Sdas	    RETURNI(x);			/* x is tiny */
67251599Sdas	SET_LDBL_EXPSIGN(x, ix);
68251599Sdas	if (ix < 0x3ffe) {		/* |x| < 0.5, or misnormal */
692116Sjkh	    t = x+x;
70251599Sdas	    t = 0.5*log1pl(t+t*x/(one-x));
71141296Sdas	} else
72251599Sdas	    t = 0.5*log1pl((x+x)/(one-x));
73251599Sdas	RETURNI((hx & 0x8000) == 0 ? t : -t);
742116Sjkh}
75