e_atanh.c revision 21673 
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 21673 1997-01-14 07:20:47Z jkh $";
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
38#ifdef __STDC__
39static const double one = 1.0, huge = 1e300;
40#else
41static double one = 1.0, huge = 1e300;
42#endif
43
44#ifdef __STDC__
45static const double zero = 0.0;
46#else
47static double zero = 0.0;
48#endif
49
50#ifdef __STDC__
51	double __ieee754_atanh(double x)
52#else
53	double __ieee754_atanh(x)
54	double x;
55#endif
56{
57	double t;
58	int32_t hx,ix;
59	u_int32_t lx;
60	EXTRACT_WORDS(hx,lx,x);
61	ix = hx&0x7fffffff;
62	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
63	    return (x-x)/(x-x);
64	if(ix==0x3ff00000)
65	    return x/zero;
66	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
67	SET_HIGH_WORD(x,ix);
68	if(ix<0x3fe00000) {		/* x < 0.5 */
69	    t = x+x;
70	    t = 0.5*log1p(t+t*x/(one-x));
71	} else
72	    t = 0.5*log1p((x+x)/(one-x));
73	if(hx>=0) return t; else return -t;
74}
75