|
|
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[] = "@(#)asincos.c 5.4 (Berkeley) 6/1/90";
27: #endif /* not lint */
28:
29: /* ASIN(X)
30: * RETURNS ARC SINE OF X
31: * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
32: * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
33: *
34: * Required system supported functions:
35: * copysign(x,y)
36: * sqrt(x)
37: *
38: * Required kernel function:
39: * atan2(y,x)
40: *
41: * Method :
42: * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
43: * computed as follows
44: * 1-x*x if x < 0.5,
45: * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
46: *
47: * Special cases:
48: * if x is NaN, return x itself;
49: * if |x|>1, return NaN.
50: *
51: * Accuracy:
52: * 1) If atan2() uses machine PI, then
53: *
54: * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
55: * and PI is the exact pi rounded to machine precision (see atan2 for
56: * details):
57: *
58: * in decimal:
59: * pi = 3.141592653589793 23846264338327 .....
60: * 53 bits PI = 3.141592653589793 115997963 ..... ,
61: * 56 bits PI = 3.141592653589793 227020265 ..... ,
62: *
63: * in hexadecimal:
64: * pi = 3.243F6A8885A308D313198A2E....
65: * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
66: * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
67: *
68: * In a test run with more than 200,000 random arguments on a VAX, the
69: * maximum observed error in ulps (units in the last place) was
70: * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x)));
71: *
72: * 2) If atan2() uses true pi, then
73: *
74: * asin(x) returns the exact asin(x) with error below about 2 ulps.
75: *
76: * In a test run with more than 1,024,000 random arguments on a VAX, the
77: * maximum observed error in ulps (units in the last place) was
78: * 1.99 ulps.
79: */
80:
81: double asin(x)
82: double x;
83: {
84: double s,t,copysign(),atan2(),sqrt(),one=1.0;
85: #if !defined(vax)&&!defined(tahoe)
86: if(x!=x) return(x); /* x is NaN */
87: #endif /* !defined(vax)&&!defined(tahoe) */
88: s=copysign(x,one);
89: if(s <= 0.5)
90: return(atan2(x,sqrt(one-x*x)));
91: else
92: { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
93:
94: }
95:
96: /* ACOS(X)
97: * RETURNS ARC COS OF X
98: * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
99: * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
100: *
101: * Required system supported functions:
102: * copysign(x,y)
103: * sqrt(x)
104: *
105: * Required kernel function:
106: * atan2(y,x)
107: *
108: * Method :
109: * ________
110: * / 1 - x
111: * acos(x) = 2*atan2( / -------- , 1 ) .
112: * \/ 1 + x
113: *
114: * Special cases:
115: * if x is NaN, return x itself;
116: * if |x|>1, return NaN.
117: *
118: * Accuracy:
119: * 1) If atan2() uses machine PI, then
120: *
121: * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
122: * and PI is the exact pi rounded to machine precision (see atan2 for
123: * details):
124: *
125: * in decimal:
126: * pi = 3.141592653589793 23846264338327 .....
127: * 53 bits PI = 3.141592653589793 115997963 ..... ,
128: * 56 bits PI = 3.141592653589793 227020265 ..... ,
129: *
130: * in hexadecimal:
131: * pi = 3.243F6A8885A308D313198A2E....
132: * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
133: * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
134: *
135: * In a test run with more than 200,000 random arguments on a VAX, the
136: * maximum observed error in ulps (units in the last place) was
137: * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x)));
138: *
139: * 2) If atan2() uses true pi, then
140: *
141: * acos(x) returns the exact acos(x) with error below about 2 ulps.
142: *
143: * In a test run with more than 1,024,000 random arguments on a VAX, the
144: * maximum observed error in ulps (units in the last place) was
145: * 2.15 ulps.
146: */
147:
148: double acos(x)
149: double x;
150: {
151: double t,copysign(),atan2(),sqrt(),one=1.0;
152: #if !defined(vax)&&!defined(tahoe)
153: if(x!=x) return(x);
154: #endif /* !defined(vax)&&!defined(tahoe) */
155: if( x != -1.0)
156: t=atan2(sqrt((one-x)/(one+x)),one);
157: else
158: t=atan2(one,0.0); /* t = PI/2 */
159: return(t+t);
160: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.