c3c2a1df04689044c5d0855e53b0fa429d8e5f26
[kopensolaris-gnu/glibc.git] / sysdeps / generic / tanh.c
1 /*
2  * Copyright (c) 1985 Regents of the University of California.
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  * 1. Redistributions of source code must retain the above copyright
9  *    notice, this list of conditions and the following disclaimer.
10  * 2. Redistributions in binary form must reproduce the above copyright
11  *    notice, this list of conditions and the following disclaimer in the
12  *    documentation and/or other materials provided with the distribution.
13  * 3. All advertising materials mentioning features or use of this software
14  *    must display the following acknowledgement:
15  *      This product includes software developed by the University of
16  *      California, Berkeley and its contributors.
17  * 4. Neither the name of the University nor the names of its contributors
18  *    may be used to endorse or promote products derived from this software
19  *    without specific prior written permission.
20  *
21  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31  * SUCH DAMAGE.
32  */
33
34 #ifndef lint
35 static char sccsid[] = "@(#)tanh.c      5.5 (Berkeley) 10/9/90";
36 #endif /* not lint */
37
38 /* TANH(X)
39  * RETURN THE HYPERBOLIC TANGENT OF X
40  * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
41  * CODED IN C BY K.C. NG, 1/8/85; 
42  * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85.
43  *
44  * Required system supported functions :
45  *      copysign(x,y)
46  *      finite(x)
47  *
48  * Required kernel function:
49  *      expm1(x)        ...exp(x)-1
50  *
51  * Method :
52  *      1. reduce x to non-negative by tanh(-x) = - tanh(x).
53  *      2.
54  *          0      <  x <=  1.e-10 :  tanh(x) := x
55  *                                                -expm1(-2x)
56  *          1.e-10 <  x <=  1      :  tanh(x) := --------------
57  *                                               expm1(-2x) + 2
58  *                                                        2
59  *          1      <= x <=  22.0   :  tanh(x) := 1 -  ---------------
60  *                                                    expm1(2x) + 2
61  *          22.0   <  x <= INF     :  tanh(x) := 1.
62  *
63  *      Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1.
64  *
65  * Special cases:
66  *      tanh(NaN) is NaN;
67  *      only tanh(0)=0 is exact for finite argument.
68  *
69  * Accuracy:
70  *      tanh(x) returns the exact hyperbolic tangent of x nealy rounded.
71  *      In a test run with 1,024,000 random arguments on a VAX, the maximum
72  *      observed error was 2.22 ulps (units in the last place).
73  */
74
75 double tanh(x)
76 double x;
77 {
78         static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10;
79         double expm1(), t, copysign(), sign;
80         int finite();
81
82 #if !defined(vax)&&!defined(tahoe)
83         if(x!=x) return(x);     /* x is NaN */
84 #endif  /* !defined(vax)&&!defined(tahoe) */
85
86         sign=copysign(one,x);
87         x=copysign(x,one);
88         if(x < 22.0) 
89             if( x > one )
90                 return(copysign(one-two/(expm1(x+x)+two),sign));
91             else if ( x > small )
92                 {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));}
93             else                /* raise the INEXACT flag for non-zero x */
94                 {big+x; return(copysign(x,sign));}
95         else if(finite(x))
96             return (sign+1.0E-37); /* raise the INEXACT flag */
97         else
98             return(sign);       /* x is +- INF */
99 }