134126Sbostic /* 224606Szliu * Copyright (c) 1985 Regents of the University of California. 334126Sbostic * All rights reserved. 434126Sbostic * 5*42657Sbostic * %sccs.include.redist.c% 634126Sbostic * 734126Sbostic * All recipients should regard themselves as participants in an ongoing 834126Sbostic * research project and hence should feel obligated to report their 934126Sbostic * experiences (good or bad) with these elementary function codes, using 1034126Sbostic * the sendbug(8) program, to the authors. 1124606Szliu */ 1224606Szliu 1324606Szliu #ifndef lint 14*42657Sbostic static char sccsid[] = "@(#)sinh.c 5.5 (Berkeley) 06/01/90"; 1534126Sbostic #endif /* not lint */ 1624606Szliu 1724606Szliu /* SINH(X) 1824606Szliu * RETURN THE HYPERBOLIC SINE OF X 1924606Szliu * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) 2024606Szliu * CODED IN C BY K.C. NG, 1/8/85; 2124606Szliu * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85. 2224606Szliu * 2324606Szliu * Required system supported functions : 2424606Szliu * copysign(x,y) 2524606Szliu * scalb(x,N) 2624606Szliu * 2724606Szliu * Required kernel functions: 2824606Szliu * expm1(x) ...return exp(x)-1 2924606Szliu * 3024606Szliu * Method : 3124606Szliu * 1. reduce x to non-negative by sinh(-x) = - sinh(x). 3224606Szliu * 2. 3324606Szliu * 3424606Szliu * expm1(x) + expm1(x)/(expm1(x)+1) 3524606Szliu * 0 <= x <= lnovfl : sinh(x) := -------------------------------- 3624606Szliu * 2 3724606Szliu * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow) 3824606Szliu * lnovfl+ln2 < x < INF : overflow to INF 3924606Szliu * 4024606Szliu * 4124606Szliu * Special cases: 4224606Szliu * sinh(x) is x if x is +INF, -INF, or NaN. 4324606Szliu * only sinh(0)=0 is exact for finite argument. 4424606Szliu * 4524606Szliu * Accuracy: 4624606Szliu * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In 4724606Szliu * a test run with 1,024,000 random arguments on a VAX, the maximum 4824606Szliu * observed error was 1.93 ulps (units in the last place). 4924606Szliu * 5024606Szliu * Constants: 5124606Szliu * The hexadecimal values are the intended ones for the following constants. 5224606Szliu * The decimal values may be used, provided that the compiler will convert 5324606Szliu * from decimal to binary accurately enough to produce the hexadecimal values 5424606Szliu * shown. 5524606Szliu */ 5624606Szliu 5735679Sbostic #include "mathimpl.h" 5835679Sbostic 5935679Sbostic vc(mln2hi, 8.8029691931113054792E1 ,0f33,43b0,2bdb,c7e2, 7, .B00F33C7E22BDB) 6035679Sbostic vc(mln2lo,-4.9650192275318476525E-16 ,1b60,a70f,582a,279e, -50,-.8F1B60279E582A) 6135679Sbostic vc(lnovfl, 8.8029691931113053016E1 ,0f33,43b0,2bda,c7e2, 7, .B00F33C7E22BDA) 6235679Sbostic 6335679Sbostic ic(mln2hi, 7.0978271289338397310E2, 10, 1.62E42FEFA39EF) 6435679Sbostic ic(mln2lo, 2.3747039373786107478E-14, -45, 1.ABC9E3B39803F) 6535679Sbostic ic(lnovfl, 7.0978271289338397310E2, 9, 1.62E42FEFA39EF) 6635679Sbostic 6735679Sbostic #ifdef vccast 6835679Sbostic #define mln2hi vccast(mln2hi) 6935679Sbostic #define mln2lo vccast(mln2lo) 7035679Sbostic #define lnovfl vccast(lnovfl) 7135679Sbostic #endif 7235679Sbostic 7331853Szliu #if defined(vax)||defined(tahoe) 7424606Szliu static max = 126 ; 7531853Szliu #else /* defined(vax)||defined(tahoe) */ 7624606Szliu static max = 1023 ; 7731853Szliu #endif /* defined(vax)||defined(tahoe) */ 7824606Szliu 7924606Szliu 8024606Szliu double sinh(x) 8124606Szliu double x; 8224606Szliu { 8335679Sbostic static const double one=1.0, half=1.0/2.0 ; 8435679Sbostic double t, sign; 8531853Szliu #if !defined(vax)&&!defined(tahoe) 8624606Szliu if(x!=x) return(x); /* x is NaN */ 8731853Szliu #endif /* !defined(vax)&&!defined(tahoe) */ 8824606Szliu sign=copysign(one,x); 8924606Szliu x=copysign(x,one); 9024606Szliu if(x<lnovfl) 9124606Szliu {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));} 9224606Szliu 9324606Szliu else if(x <= lnovfl+0.7) 9424606Szliu /* subtract x by ln(2^(max+1)) and return 2^max*exp(x) 9524606Szliu to avoid unnecessary overflow */ 9624606Szliu return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign)); 9724606Szliu 9824606Szliu else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */ 9924606Szliu return( expm1(x)*sign ); 10024606Szliu } 101