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%
624606Szliu */
724606Szliu
824606Szliu #ifndef lint
9*61309Sbostic static char sccsid[] = "@(#)sinh.c 8.1 (Berkeley) 06/04/93";
1034126Sbostic #endif /* not lint */
1124606Szliu
1224606Szliu /* SINH(X)
1324606Szliu * RETURN THE HYPERBOLIC SINE OF X
1424606Szliu * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
1524606Szliu * CODED IN C BY K.C. NG, 1/8/85;
1624606Szliu * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85.
1724606Szliu *
1824606Szliu * Required system supported functions :
1924606Szliu * copysign(x,y)
2024606Szliu * scalb(x,N)
2124606Szliu *
2224606Szliu * Required kernel functions:
2324606Szliu * expm1(x) ...return exp(x)-1
2424606Szliu *
2524606Szliu * Method :
2624606Szliu * 1. reduce x to non-negative by sinh(-x) = - sinh(x).
2724606Szliu * 2.
2824606Szliu *
2924606Szliu * expm1(x) + expm1(x)/(expm1(x)+1)
3024606Szliu * 0 <= x <= lnovfl : sinh(x) := --------------------------------
3124606Szliu * 2
3224606Szliu * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow)
3324606Szliu * lnovfl+ln2 < x < INF : overflow to INF
3424606Szliu *
3524606Szliu *
3624606Szliu * Special cases:
3724606Szliu * sinh(x) is x if x is +INF, -INF, or NaN.
3824606Szliu * only sinh(0)=0 is exact for finite argument.
3924606Szliu *
4024606Szliu * Accuracy:
4124606Szliu * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In
4224606Szliu * a test run with 1,024,000 random arguments on a VAX, the maximum
4324606Szliu * observed error was 1.93 ulps (units in the last place).
4424606Szliu *
4524606Szliu * Constants:
4624606Szliu * The hexadecimal values are the intended ones for the following constants.
4724606Szliu * The decimal values may be used, provided that the compiler will convert
4824606Szliu * from decimal to binary accurately enough to produce the hexadecimal values
4924606Szliu * shown.
5024606Szliu */
5124606Szliu
5235679Sbostic #include "mathimpl.h"
5335679Sbostic
5435679Sbostic vc(mln2hi, 8.8029691931113054792E1 ,0f33,43b0,2bdb,c7e2, 7, .B00F33C7E22BDB)
5535679Sbostic vc(mln2lo,-4.9650192275318476525E-16 ,1b60,a70f,582a,279e, -50,-.8F1B60279E582A)
5635679Sbostic vc(lnovfl, 8.8029691931113053016E1 ,0f33,43b0,2bda,c7e2, 7, .B00F33C7E22BDA)
5735679Sbostic
5835679Sbostic ic(mln2hi, 7.0978271289338397310E2, 10, 1.62E42FEFA39EF)
5935679Sbostic ic(mln2lo, 2.3747039373786107478E-14, -45, 1.ABC9E3B39803F)
6035679Sbostic ic(lnovfl, 7.0978271289338397310E2, 9, 1.62E42FEFA39EF)
6135679Sbostic
6235679Sbostic #ifdef vccast
6335679Sbostic #define mln2hi vccast(mln2hi)
6435679Sbostic #define mln2lo vccast(mln2lo)
6535679Sbostic #define lnovfl vccast(lnovfl)
6635679Sbostic #endif
6735679Sbostic
6831853Szliu #if defined(vax)||defined(tahoe)
6924606Szliu static max = 126 ;
7031853Szliu #else /* defined(vax)||defined(tahoe) */
7124606Szliu static max = 1023 ;
7231853Szliu #endif /* defined(vax)||defined(tahoe) */
7324606Szliu
7424606Szliu
sinh(x)7524606Szliu double sinh(x)
7624606Szliu double x;
7724606Szliu {
7835679Sbostic static const double one=1.0, half=1.0/2.0 ;
7935679Sbostic double t, sign;
8031853Szliu #if !defined(vax)&&!defined(tahoe)
8124606Szliu if(x!=x) return(x); /* x is NaN */
8231853Szliu #endif /* !defined(vax)&&!defined(tahoe) */
8324606Szliu sign=copysign(one,x);
8424606Szliu x=copysign(x,one);
8524606Szliu if(x<lnovfl)
8624606Szliu {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));}
8724606Szliu
8824606Szliu else if(x <= lnovfl+0.7)
8924606Szliu /* subtract x by ln(2^(max+1)) and return 2^max*exp(x)
9024606Szliu to avoid unnecessary overflow */
9124606Szliu return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign));
9224606Szliu
9324606Szliu else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */
9424606Szliu return( expm1(x)*sign );
9524606Szliu }
96