|
|
1.1 root 1: /*
2: * Copyright (c) 1985 Regents of the University of California.
3: *
4: * Use and reproduction of this software are granted in accordance with
5: * the terms and conditions specified in the Berkeley Software License
6: * Agreement (in particular, this entails acknowledgement of the programs'
7: * source, and inclusion of this notice) with the additional understanding
8: * that all recipients should regard themselves as participants in an
9: * ongoing research project and hence should feel obligated to report
10: * their experiences (good or bad) with these elementary function codes,
11: * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
12: */
13:
14: #ifndef lint
15: static char sccsid[] = "@(#)asincos.c 1.1 (Berkeley) 8/21/85";
16: #endif not lint
17:
18: /* ASIN(X)
19: * RETURNS ARC SINE OF X
20: * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
21: * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
22: *
23: * Required system supported functions:
24: * copysign(x,y)
25: * sqrt(x)
26: *
27: * Required kernel function:
28: * atan2(y,x)
29: *
30: * Method :
31: * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
32: * computed as follows
33: * 1-x*x if x < 0.5,
34: * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
35: *
36: * Special cases:
37: * if x is NaN, return x itself;
38: * if |x|>1, return NaN.
39: *
40: * Accuracy:
41: * 1) If atan2() uses machine PI, then
42: *
43: * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
44: * and PI is the exact pi rounded to machine precision (see atan2 for
45: * details):
46: *
47: * in decimal:
48: * pi = 3.141592653589793 23846264338327 .....
49: * 53 bits PI = 3.141592653589793 115997963 ..... ,
50: * 56 bits PI = 3.141592653589793 227020265 ..... ,
51: *
52: * in hexadecimal:
53: * pi = 3.243F6A8885A308D313198A2E....
54: * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
55: * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
56: *
57: * In a test run with more than 200,000 random arguments on a VAX, the
58: * maximum observed error in ulps (units in the last place) was
59: * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x)));
60: *
61: * 2) If atan2() uses true pi, then
62: *
63: * asin(x) returns the exact asin(x) with error below about 2 ulps.
64: *
65: * In a test run with more than 1,024,000 random arguments on a VAX, the
66: * maximum observed error in ulps (units in the last place) was
67: * 1.99 ulps.
68: */
69:
70: double asin(x)
71: double x;
72: {
73: double s,t,copysign(),atan2(),sqrt(),one=1.0;
74: #ifndef VAX
75: if(x!=x) return(x); /* x is NaN */
76: #endif
77: s=copysign(x,one);
78: if(s <= 0.5)
79: return(atan2(x,sqrt(one-x*x)));
80: else
81: { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
82:
83: }
84:
85: /* ACOS(X)
86: * RETURNS ARC COS OF X
87: * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
88: * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
89: *
90: * Required system supported functions:
91: * copysign(x,y)
92: * sqrt(x)
93: *
94: * Required kernel function:
95: * atan2(y,x)
96: *
97: * Method :
98: * ________
99: * / 1 - x
100: * acos(x) = 2*atan2( / -------- , 1 ) .
101: * \/ 1 + x
102: *
103: * Special cases:
104: * if x is NaN, return x itself;
105: * if |x|>1, return NaN.
106: *
107: * Accuracy:
108: * 1) If atan2() uses machine PI, then
109: *
110: * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
111: * and PI is the exact pi rounded to machine precision (see atan2 for
112: * details):
113: *
114: * in decimal:
115: * pi = 3.141592653589793 23846264338327 .....
116: * 53 bits PI = 3.141592653589793 115997963 ..... ,
117: * 56 bits PI = 3.141592653589793 227020265 ..... ,
118: *
119: * in hexadecimal:
120: * pi = 3.243F6A8885A308D313198A2E....
121: * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
122: * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
123: *
124: * In a test run with more than 200,000 random arguments on a VAX, the
125: * maximum observed error in ulps (units in the last place) was
126: * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x)));
127: *
128: * 2) If atan2() uses true pi, then
129: *
130: * acos(x) returns the exact acos(x) with error below about 2 ulps.
131: *
132: * In a test run with more than 1,024,000 random arguments on a VAX, the
133: * maximum observed error in ulps (units in the last place) was
134: * 2.15 ulps.
135: */
136:
137: double acos(x)
138: double x;
139: {
140: double t,copysign(),atan2(),sqrt(),one=1.0;
141: #ifndef VAX
142: if(x!=x) return(x);
143: #endif
144: if( x != -1.0)
145: t=atan2(sqrt((one-x)/(one+x)),one);
146: else
147: t=atan2(one,0.0); /* t = PI/2 */
148: return(t+t);
149: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.