1 /* 2 * Copyright (c) 1985 Regents of the University of California. 3 * 4 * Use and reproduction of this software are granted in accordance with 5 * the terms and conditions specified in the Berkeley Software License 6 * Agreement (in particular, this entails acknowledgement of the programs' 7 * source, and inclusion of this notice) with the additional understanding 8 * that all recipients should regard themselves as participants in an 9 * ongoing research project and hence should feel obligated to report 10 * their experiences (good or bad) with these elementary function codes, 11 * using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12 */ 13 14 #ifndef lint 15 static char sccsid[] = 16 "@(#)atanh.c 1.2 (Berkeley) 8/21/85; 5.1 (ucb.elefunt) 11/30/87"; 17 #endif /* not lint */ 18 19 /* ATANH(X) 20 * RETURN THE HYPERBOLIC ARC TANGENT OF X 21 * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 22 * CODED IN C BY K.C. NG, 1/8/85; 23 * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85. 24 * 25 * Required kernel function: 26 * log1p(x) ...return log(1+x) 27 * 28 * Method : 29 * Return 30 * 1 2x x 31 * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) 32 * 2 1 - x 1 - x 33 * 34 * Special cases: 35 * atanh(x) is NaN if |x| > 1 with signal; 36 * atanh(NaN) is that NaN with no signal; 37 * atanh(+-1) is +-INF with signal. 38 * 39 * Accuracy: 40 * atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded. 41 * In a test run with 512,000 random arguments on a VAX, the maximum 42 * observed error was 1.87 ulps (units in the last place) at 43 * x= -3.8962076028810414000e-03. 44 */ 45 #if defined(vax)||defined(tahoe) 46 #include <errno.h> 47 #endif /* defined(vax)||defined(tahoe) */ 48 49 double atanh(x) 50 double x; 51 { 52 double copysign(),log1p(),z; 53 z = copysign(0.5,x); 54 x = copysign(x,1.0); 55 #if defined(vax)||defined(tahoe) 56 if (x == 1.0) { 57 extern double infnan(); 58 return(copysign(1.0,z)*infnan(ERANGE)); /* sign(x)*INF */ 59 } 60 #endif /* defined(vax)||defined(tahoe) */ 61 x = x/(1.0-x); 62 return( z*log1p(x+x) ); 63 } 64