|
|
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[] = "@(#)asinh.c 5.5 (Berkeley) 6/1/90";
27: #endif /* not lint */
28:
29: /* ASINH(X)
30: * RETURN THE INVERSE HYPERBOLIC SINE OF X
31: * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
32: * CODED IN C BY K.C. NG, 2/16/85;
33: * REVISED BY K.C. NG on 3/7/85, 3/24/85, 4/16/85.
34: *
35: * Required system supported functions :
36: * copysign(x,y)
37: * sqrt(x)
38: *
39: * Required kernel function:
40: * log1p(x) ...return log(1+x)
41: *
42: * Method :
43: * Based on
44: * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
45: * we have
46: * asinh(x) := x if 1+x*x=1,
47: * := sign(x)*(log1p(x)+ln2)) if sqrt(1+x*x)=x, else
48: * := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) )
49: *
50: * Accuracy:
51: * asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded.
52: * In a test run with 52,000 random arguments on a VAX, the maximum
53: * observed error was 1.58 ulps (units in the last place).
54: *
55: * Constants:
56: * The hexadecimal values are the intended ones for the following constants.
57: * The decimal values may be used, provided that the compiler will convert
58: * from decimal to binary accurately enough to produce the hexadecimal values
59: * shown.
60: */
61: #include "mathimpl.h"
62:
63: vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000)
64: vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
65:
66: ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000)
67: ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76)
68:
69: #ifdef vccast
70: #define ln2hi vccast(ln2hi)
71: #define ln2lo vccast(ln2lo)
72: #endif
73:
74: double asinh(x)
75: double x;
76: {
77: double t,s;
78: const static double small=1.0E-10, /* fl(1+small*small) == 1 */
79: big =1.0E20, /* fl(1+big) == big */
80: one =1.0 ;
81:
82: #if !defined(vax)&&!defined(tahoe)
83: if(x!=x) return(x); /* x is NaN */
84: #endif /* !defined(vax)&&!defined(tahoe) */
85: if((t=copysign(x,one))>small)
86: if(t<big) {
87: s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); }
88: else /* if |x| > big */
89: {s=log1p(t)+ln2lo; return(copysign(s+ln2hi,x));}
90: else /* if |x| < small */
91: return(x);
92: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.