|
|
1.1 root 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 are permitted
6: * provided that: (1) source distributions retain this entire copyright
7: * notice and comment, and (2) distributions including binaries display
8: * the following acknowledgement: ``This product includes software
9: * developed by the University of California, Berkeley and its contributors''
10: * in the documentation or other materials provided with the distribution
11: * and in all advertising materials mentioning features or use of this
12: * software. Neither the name of the University nor the names of its
13: * contributors may be used to endorse or promote products derived
14: * from this software without specific prior written permission.
15: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
16: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
17: * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
18: *
19: * All recipients should regard themselves as participants in an ongoing
20: * research project and hence should feel obligated to report their
21: * experiences (good or bad) with these elementary function codes, using
22: * the sendbug(8) program, to the authors.
23: */
24:
25: #ifndef lint
26: static char sccsid[] = "@(#)tanh.c 5.4 (Berkeley) 6/1/90";
27: #endif /* not lint */
28:
29: /* TANH(X)
30: * RETURN THE HYPERBOLIC TANGENT OF X
31: * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
32: * CODED IN C BY K.C. NG, 1/8/85;
33: * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85.
34: *
35: * Required system supported functions :
36: * copysign(x,y)
37: * finite(x)
38: *
39: * Required kernel function:
40: * expm1(x) ...exp(x)-1
41: *
42: * Method :
43: * 1. reduce x to non-negative by tanh(-x) = - tanh(x).
44: * 2.
45: * 0 < x <= 1.e-10 : tanh(x) := x
46: * -expm1(-2x)
47: * 1.e-10 < x <= 1 : tanh(x) := --------------
48: * expm1(-2x) + 2
49: * 2
50: * 1 <= x <= 22.0 : tanh(x) := 1 - ---------------
51: * expm1(2x) + 2
52: * 22.0 < x <= INF : tanh(x) := 1.
53: *
54: * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1.
55: *
56: * Special cases:
57: * tanh(NaN) is NaN;
58: * only tanh(0)=0 is exact for finite argument.
59: *
60: * Accuracy:
61: * tanh(x) returns the exact hyperbolic tangent of x nealy rounded.
62: * In a test run with 1,024,000 random arguments on a VAX, the maximum
63: * observed error was 2.22 ulps (units in the last place).
64: */
65:
66: double tanh(x)
67: double x;
68: {
69: static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10;
70: double expm1(), t, copysign(), sign;
71: int finite();
72:
73: #if !defined(vax)&&!defined(tahoe)
74: if(x!=x) return(x); /* x is NaN */
75: #endif /* !defined(vax)&&!defined(tahoe) */
76:
77: sign=copysign(one,x);
78: x=copysign(x,one);
79: if(x < 22.0)
80: if( x > one )
81: return(copysign(one-two/(expm1(x+x)+two),sign));
82: else if ( x > small )
83: {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));}
84: else /* raise the INEXACT flag for non-zero x */
85: {big+x; return(copysign(x,sign));}
86: else if(finite(x))
87: return (sign+1.0E-37); /* raise the INEXACT flag */
88: else
89: return(sign); /* x is +- INF */
90: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.