|
|
1.1 ! root 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[] = "@(#)tanh.c 4.3 (Berkeley) 8/21/85"; ! 16: #endif not lint ! 17: ! 18: /* TANH(X) ! 19: * RETURN THE HYPERBOLIC TANGENT OF X ! 20: * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) ! 21: * CODED IN C BY K.C. NG, 1/8/85; ! 22: * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85. ! 23: * ! 24: * Required system supported functions : ! 25: * copysign(x,y) ! 26: * finite(x) ! 27: * ! 28: * Required kernel function: ! 29: * expm1(x) ...exp(x)-1 ! 30: * ! 31: * Method : ! 32: * 1. reduce x to non-negative by tanh(-x) = - tanh(x). ! 33: * 2. ! 34: * 0 < x <= 1.e-10 : tanh(x) := x ! 35: * -expm1(-2x) ! 36: * 1.e-10 < x <= 1 : tanh(x) := -------------- ! 37: * expm1(-2x) + 2 ! 38: * 2 ! 39: * 1 <= x <= 22.0 : tanh(x) := 1 - --------------- ! 40: * expm1(2x) + 2 ! 41: * 22.0 < x <= INF : tanh(x) := 1. ! 42: * ! 43: * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1. ! 44: * ! 45: * Special cases: ! 46: * tanh(NaN) is NaN; ! 47: * only tanh(0)=0 is exact for finite argument. ! 48: * ! 49: * Accuracy: ! 50: * tanh(x) returns the exact hyperbolic tangent of x nealy rounded. ! 51: * In a test run with 1,024,000 random arguments on a VAX, the maximum ! 52: * observed error was 2.22 ulps (units in the last place). ! 53: */ ! 54: ! 55: double tanh(x) ! 56: double x; ! 57: { ! 58: static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10; ! 59: double expm1(), t, copysign(), sign; ! 60: int finite(); ! 61: ! 62: #ifndef VAX ! 63: if(x!=x) return(x); /* x is NaN */ ! 64: #endif ! 65: ! 66: sign=copysign(one,x); ! 67: x=copysign(x,one); ! 68: if(x < 22.0) ! 69: if( x > one ) ! 70: return(copysign(one-two/(expm1(x+x)+two),sign)); ! 71: else if ( x > small ) ! 72: {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));} ! 73: else /* raise the INEXACT flag for non-zero x */ ! 74: {big+x; return(copysign(x,sign));} ! 75: else if(finite(x)) ! 76: return (sign+1.0E-37); /* raise the INEXACT flag */ ! 77: else ! 78: return(sign); /* x is +- INF */ ! 79: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.