134124Sbostic /*
2*61285Sbostic * Copyright (c) 1985, 1993
3*61285Sbostic * The Regents of the University of California. All rights reserved.
434124Sbostic *
542657Sbostic * %sccs.include.redist.c%
624590Szliu */
724590Szliu
824590Szliu #ifndef lint
9*61285Sbostic static char sccsid[] = "@(#)atanh.c 8.1 (Berkeley) 06/04/93";
1034124Sbostic #endif /* not lint */
1124590Szliu
1224590Szliu /* ATANH(X)
1324590Szliu * RETURN THE HYPERBOLIC ARC TANGENT OF X
1424590Szliu * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
1524590Szliu * CODED IN C BY K.C. NG, 1/8/85;
1624590Szliu * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85.
1724590Szliu *
1824590Szliu * Required kernel function:
1924590Szliu * log1p(x) ...return log(1+x)
2024590Szliu *
2124590Szliu * Method :
2224590Szliu * Return
2324590Szliu * 1 2x x
2424590Szliu * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
2524590Szliu * 2 1 - x 1 - x
2624590Szliu *
2724590Szliu * Special cases:
2824590Szliu * atanh(x) is NaN if |x| > 1 with signal;
2924590Szliu * atanh(NaN) is that NaN with no signal;
3024590Szliu * atanh(+-1) is +-INF with signal.
3124590Szliu *
3224590Szliu * Accuracy:
3324590Szliu * atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded.
3424590Szliu * In a test run with 512,000 random arguments on a VAX, the maximum
3524590Szliu * observed error was 1.87 ulps (units in the last place) at
3624590Szliu * x= -3.8962076028810414000e-03.
3724590Szliu */
3835679Sbostic #include "mathimpl.h"
3935679Sbostic
4031853Szliu #if defined(vax)||defined(tahoe)
4124590Szliu #include <errno.h>
4231853Szliu #endif /* defined(vax)||defined(tahoe) */
4324590Szliu
atanh(x)4424590Szliu double atanh(x)
4524590Szliu double x;
4624590Szliu {
4735679Sbostic double z;
4824590Szliu z = copysign(0.5,x);
4924590Szliu x = copysign(x,1.0);
5031853Szliu #if defined(vax)||defined(tahoe)
5124590Szliu if (x == 1.0) {
5224590Szliu return(copysign(1.0,z)*infnan(ERANGE)); /* sign(x)*INF */
5324590Szliu }
5431853Szliu #endif /* defined(vax)||defined(tahoe) */
5524590Szliu x = x/(1.0-x);
5624590Szliu return( z*log1p(x+x) );
5724590Szliu }
58