e_atanhl.c revision 97407 
1/* @(#)e_atanh.c 5.1 93/09/24 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#ifndef lint
14static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_atanh.c 97407 2002-05-28 17:03:12Z alfred $";
15#endif
16
17/* __ieee754_atanh(x)
18 * Method :
19 *    1.Reduced x to positive by atanh(-x) = -atanh(x)
20 *    2.For x>=0.5
21 *                  1              2x                          x
22 *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
23 *                  2             1 - x                      1 - x
24 *
25 * 	For x<0.5
26 *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
27 *
28 * Special cases:
29 *	atanh(x) is NaN if |x| > 1 with signal;
30 *	atanh(NaN) is that NaN with no signal;
31 *	atanh(+-1) is +-INF with signal.
32 *
33 */
34
35#include "math.h"
36#include "math_private.h"
37
38static const double one = 1.0, huge = 1e300;
39static const double zero = 0.0;
40
41double
42__ieee754_atanh(double x)
43{
44	double t;
45	int32_t hx,ix;
46	u_int32_t lx;
47	EXTRACT_WORDS(hx,lx,x);
48	ix = hx&0x7fffffff;
49	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
50	    return (x-x)/(x-x);
51	if(ix==0x3ff00000)
52	    return x/zero;
53	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
54	SET_HIGH_WORD(x,ix);
55	if(ix<0x3fe00000) {		/* x < 0.5 */
56	    t = x+x;
57	    t = 0.5*log1p(t+t*x/(one-x));
58	} else
59	    t = 0.5*log1p((x+x)/(one-x));
60	if(hx>=0) return t; else return -t;
61}
62