134126Sbostic /*
2*61309Sbostic * Copyright (c) 1985, 1993
3*61309Sbostic * The Regents of the University of California. All rights reserved.
434126Sbostic *
542657Sbostic * %sccs.include.redist.c%
624607Szliu */
724607Szliu
824607Szliu #ifndef lint
9*61309Sbostic static char sccsid[] = "@(#)tanh.c 8.1 (Berkeley) 06/04/93";
1034126Sbostic #endif /* not lint */
1124607Szliu
1224607Szliu /* TANH(X)
1324607Szliu * RETURN THE HYPERBOLIC TANGENT OF X
1424607Szliu * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
1524607Szliu * CODED IN C BY K.C. NG, 1/8/85;
1624607Szliu * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85.
1724607Szliu *
1824607Szliu * Required system supported functions :
1924607Szliu * copysign(x,y)
2024607Szliu * finite(x)
2124607Szliu *
2224607Szliu * Required kernel function:
2324607Szliu * expm1(x) ...exp(x)-1
2424607Szliu *
2524607Szliu * Method :
2624607Szliu * 1. reduce x to non-negative by tanh(-x) = - tanh(x).
2724607Szliu * 2.
2824607Szliu * 0 < x <= 1.e-10 : tanh(x) := x
2924607Szliu * -expm1(-2x)
3024607Szliu * 1.e-10 < x <= 1 : tanh(x) := --------------
3124607Szliu * expm1(-2x) + 2
3224607Szliu * 2
3324607Szliu * 1 <= x <= 22.0 : tanh(x) := 1 - ---------------
3424607Szliu * expm1(2x) + 2
3524607Szliu * 22.0 < x <= INF : tanh(x) := 1.
3624607Szliu *
3724607Szliu * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1.
3824607Szliu *
3924607Szliu * Special cases:
4024607Szliu * tanh(NaN) is NaN;
4124607Szliu * only tanh(0)=0 is exact for finite argument.
4224607Szliu *
4324607Szliu * Accuracy:
4424607Szliu * tanh(x) returns the exact hyperbolic tangent of x nealy rounded.
4524607Szliu * In a test run with 1,024,000 random arguments on a VAX, the maximum
4624607Szliu * observed error was 2.22 ulps (units in the last place).
4724607Szliu */
4824607Szliu
tanh(x)4924607Szliu double tanh(x)
5024607Szliu double x;
5124607Szliu {
5224607Szliu static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10;
5324607Szliu double expm1(), t, copysign(), sign;
5424607Szliu int finite();
5524607Szliu
5631853Szliu #if !defined(vax)&&!defined(tahoe)
5724607Szliu if(x!=x) return(x); /* x is NaN */
5831853Szliu #endif /* !defined(vax)&&!defined(tahoe) */
5924607Szliu
6024607Szliu sign=copysign(one,x);
6124607Szliu x=copysign(x,one);
6224607Szliu if(x < 22.0)
6324607Szliu if( x > one )
6424607Szliu return(copysign(one-two/(expm1(x+x)+two),sign));
6524607Szliu else if ( x > small )
6624607Szliu {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));}
6724607Szliu else /* raise the INEXACT flag for non-zero x */
6824607Szliu {big+x; return(copysign(x,sign));}
6924607Szliu else if(finite(x))
7024607Szliu return (sign+1.0E-37); /* raise the INEXACT flag */
7124607Szliu else
7224607Szliu return(sign); /* x is +- INF */
7324607Szliu }
74