Commit | Line | Data |
---|---|---|
4c182be8 | 1 | /* |
3cd1de16 | 2 | * Copyright (c) 1985 Regents of the University of California. |
4c182be8 KB |
3 | * All rights reserved. |
4 | * | |
5 | * Redistribution and use in source and binary forms are permitted | |
a399f6c8 KB |
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. | |
4c182be8 KB |
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. | |
3cd1de16 ZAL |
21 | */ |
22 | ||
23 | #ifndef lint | |
a399f6c8 | 24 | static char sccsid[] = "@(#)asincos.c 5.3 (Berkeley) %G%"; |
4c182be8 | 25 | #endif /* not lint */ |
3cd1de16 ZAL |
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; | |
859dc438 | 83 | #if !defined(vax)&&!defined(tahoe) |
3cd1de16 | 84 | if(x!=x) return(x); /* x is NaN */ |
859dc438 | 85 | #endif /* !defined(vax)&&!defined(tahoe) */ |
3cd1de16 ZAL |
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; | |
859dc438 | 150 | #if !defined(vax)&&!defined(tahoe) |
3cd1de16 | 151 | if(x!=x) return(x); |
859dc438 | 152 | #endif /* !defined(vax)&&!defined(tahoe) */ |
3cd1de16 ZAL |
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 | } |