121035Smiriam /*
221035Smiriam * Copyright (c) 1985 Regents of the University of California.
321035Smiriam *
421035Smiriam * Use and reproduction of this software are granted in accordance with
521035Smiriam * the terms and conditions specified in the Berkeley Software License
621035Smiriam * Agreement (in particular, this entails acknowledgement of the programs'
721035Smiriam * source, and inclusion of this notice) with the additional understanding
821035Smiriam * that all recipients should regard themselves as participants in an
921035Smiriam * ongoing research project and hence should feel obligated to report
1021035Smiriam * their experiences (good or bad) with these elementary function codes,
1121035Smiriam * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
1221035Smiriam */
139942Ssam
1421035Smiriam #ifndef lint
15*24349Smiriam static char sccsid[] = "@(#)tanh.c 4.3 (Berkeley) 08/21/85";
1621035Smiriam #endif not lint
179942Ssam
1821035Smiriam /* TANH(X)
1921035Smiriam * RETURN THE HYPERBOLIC TANGENT OF X
2021035Smiriam * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
2121035Smiriam * CODED IN C BY K.C. NG, 1/8/85;
2221035Smiriam * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85.
2321035Smiriam *
2421035Smiriam * Required system supported functions :
2521035Smiriam * copysign(x,y)
2621035Smiriam * finite(x)
2721035Smiriam *
2821035Smiriam * Required kernel function:
29*24349Smiriam * expm1(x) ...exp(x)-1
3021035Smiriam *
3121035Smiriam * Method :
3221035Smiriam * 1. reduce x to non-negative by tanh(-x) = - tanh(x).
33*24349Smiriam * 2.
34*24349Smiriam * 0 < x <= 1.e-10 : tanh(x) := x
35*24349Smiriam * -expm1(-2x)
36*24349Smiriam * 1.e-10 < x <= 1 : tanh(x) := --------------
37*24349Smiriam * expm1(-2x) + 2
38*24349Smiriam * 2
39*24349Smiriam * 1 <= x <= 22.0 : tanh(x) := 1 - ---------------
40*24349Smiriam * expm1(2x) + 2
41*24349Smiriam * 22.0 < x <= INF : tanh(x) := 1.
4221035Smiriam *
43*24349Smiriam * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1.
4421035Smiriam *
4521035Smiriam * Special cases:
46*24349Smiriam * tanh(NaN) is NaN;
4721035Smiriam * only tanh(0)=0 is exact for finite argument.
4821035Smiriam *
4921035Smiriam * Accuracy:
5021035Smiriam * tanh(x) returns the exact hyperbolic tangent of x nealy rounded.
5121035Smiriam * In a test run with 1,024,000 random arguments on a VAX, the maximum
5221035Smiriam * observed error was 2.22 ulps (units in the last place).
5321035Smiriam */
549942Ssam
tanh(x)5521035Smiriam double tanh(x)
5621035Smiriam double x;
579942Ssam {
58*24349Smiriam static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10;
59*24349Smiriam double expm1(), t, copysign(), sign;
6021035Smiriam int finite();
619942Ssam
62*24349Smiriam #ifndef VAX
63*24349Smiriam if(x!=x) return(x); /* x is NaN */
64*24349Smiriam #endif
659942Ssam
6621035Smiriam sign=copysign(one,x);
6721035Smiriam x=copysign(x,one);
6821035Smiriam if(x < 22.0)
6921035Smiriam if( x > one )
70*24349Smiriam return(copysign(one-two/(expm1(x+x)+two),sign));
71*24349Smiriam else if ( x > small )
72*24349Smiriam {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));}
73*24349Smiriam else /* raise the INEXACT flag for non-zero x */
74*24349Smiriam {big+x; return(copysign(x,sign));}
7521035Smiriam else if(finite(x))
76*24349Smiriam return (sign+1.0E-37); /* raise the INEXACT flag */
7721035Smiriam else
7821035Smiriam return(sign); /* x is +- INF */
799942Ssam }
80