|
|
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 the above copyright notice and this paragraph are
7: * duplicated in all such forms and that any documentation,
8: * advertising materials, and other materials related to such
9: * distribution and use acknowledge that the software was developed
10: * by the University of California, Berkeley. The name of the
11: * University may not be used to endorse or promote products derived
12: * from this software without specific prior written permission.
13: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
14: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
15: * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16: *
17: * All recipients should regard themselves as participants in an ongoing
18: * research project and hence should feel obligated to report their
19: * experiences (good or bad) with these elementary function codes, using
20: * the sendbug(8) program, to the authors.
21: */
22:
23: #ifndef lint
24: static char sccsid[] = "@(#)sinh.c 5.3 (Berkeley) 6/30/88";
25: #endif /* not lint */
26:
27: /* SINH(X)
28: * RETURN THE HYPERBOLIC SINE OF X
29: * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
30: * CODED IN C BY K.C. NG, 1/8/85;
31: * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85.
32: *
33: * Required system supported functions :
34: * copysign(x,y)
35: * scalb(x,N)
36: *
37: * Required kernel functions:
38: * expm1(x) ...return exp(x)-1
39: *
40: * Method :
41: * 1. reduce x to non-negative by sinh(-x) = - sinh(x).
42: * 2.
43: *
44: * expm1(x) + expm1(x)/(expm1(x)+1)
45: * 0 <= x <= lnovfl : sinh(x) := --------------------------------
46: * 2
47: * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow)
48: * lnovfl+ln2 < x < INF : overflow to INF
49: *
50: *
51: * Special cases:
52: * sinh(x) is x if x is +INF, -INF, or NaN.
53: * only sinh(0)=0 is exact for finite argument.
54: *
55: * Accuracy:
56: * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In
57: * a test run with 1,024,000 random arguments on a VAX, the maximum
58: * observed error was 1.93 ulps (units in the last place).
59: *
60: * Constants:
61: * The hexadecimal values are the intended ones for the following constants.
62: * The decimal values may be used, provided that the compiler will convert
63: * from decimal to binary accurately enough to produce the hexadecimal values
64: * shown.
65: */
66: #if defined(vax)||defined(tahoe)
67: #ifdef vax
68: #define _0x(A,B) 0x/**/A/**/B
69: #else /* vax */
70: #define _0x(A,B) 0x/**/B/**/A
71: #endif /* vax */
72: /* static double */
73: /* mln2hi = 8.8029691931113054792E1 , Hex 2^ 7 * .B00F33C7E22BDB */
74: /* mln2lo = -4.9650192275318476525E-16 , Hex 2^-50 * -.8F1B60279E582A */
75: /* lnovfl = 8.8029691931113053016E1 ; Hex 2^ 7 * .B00F33C7E22BDA */
76: static long mln2hix[] = { _0x(0f33,43b0), _0x(2bdb,c7e2)};
77: static long mln2lox[] = { _0x(1b60,a70f), _0x(582a,279e)};
78: static long lnovflx[] = { _0x(0f33,43b0), _0x(2bda,c7e2)};
79: #define mln2hi (*(double*)mln2hix)
80: #define mln2lo (*(double*)mln2lox)
81: #define lnovfl (*(double*)lnovflx)
82: #else /* defined(vax)||defined(tahoe) */
83: static double
84: mln2hi = 7.0978271289338397310E2 , /*Hex 2^ 10 * 1.62E42FEFA39EF */
85: mln2lo = 2.3747039373786107478E-14 , /*Hex 2^-45 * 1.ABC9E3B39803F */
86: lnovfl = 7.0978271289338397310E2 ; /*Hex 2^ 9 * 1.62E42FEFA39EF */
87: #endif /* defined(vax)||defined(tahoe) */
88:
89: #if defined(vax)||defined(tahoe)
90: static max = 126 ;
91: #else /* defined(vax)||defined(tahoe) */
92: static max = 1023 ;
93: #endif /* defined(vax)||defined(tahoe) */
94:
95:
96: double sinh(x)
97: double x;
98: {
99: static double one=1.0, half=1.0/2.0 ;
100: double expm1(), t, scalb(), copysign(), sign;
101: #if !defined(vax)&&!defined(tahoe)
102: if(x!=x) return(x); /* x is NaN */
103: #endif /* !defined(vax)&&!defined(tahoe) */
104: sign=copysign(one,x);
105: x=copysign(x,one);
106: if(x<lnovfl)
107: {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));}
108:
109: else if(x <= lnovfl+0.7)
110: /* subtract x by ln(2^(max+1)) and return 2^max*exp(x)
111: to avoid unnecessary overflow */
112: return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign));
113:
114: else /* sinh(+-INF) = +-INF, sinh(+-big no.) overflow to +-INF */
115: return( expm1(x)*sign );
116: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.