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