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, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 * 3. All advertising materials mentioning features or use of this software
14 * must display the following acknowledgement:
15 * This product includes software developed by the University of
16 * California, Berkeley and its contributors.
17 * 4. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31 * SUCH DAMAGE.
32 */
33
34 #ifndef lint
35 static char sccsid[] = "@(#)asincos.c 5.5 (Berkeley) 10/9/90";
36 #endif /* not lint */
37
38 /* ASIN(X)
39 * RETURNS ARC SINE OF X
40 * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
41 * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
42 *
43 * Required system supported functions:
44 * copysign(x,y)
45 * sqrt(x)
46 *
47 * Required kernel function:
48 * atan2(y,x)
49 *
50 * Method :
51 * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
52 * computed as follows
53 * 1-x*x if x < 0.5,
54 * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
55 *
56 * Special cases:
57 * if x is NaN, return x itself;
58 * if |x|>1, return NaN.
59 *
60 * Accuracy:
61 * 1) If atan2() uses machine PI, then
62 *
63 * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
64 * and PI is the exact pi rounded to machine precision (see atan2 for
65 * details):
66 *
67 * in decimal:
68 * pi = 3.141592653589793 23846264338327 .....
69 * 53 bits PI = 3.141592653589793 115997963 ..... ,
70 * 56 bits PI = 3.141592653589793 227020265 ..... ,
71 *
72 * in hexadecimal:
73 * pi = 3.243F6A8885A308D313198A2E....
74 * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
75 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
76 *
77 * In a test run with more than 200,000 random arguments on a VAX, the
78 * maximum observed error in ulps (units in the last place) was
79 * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x)));
80 *
81 * 2) If atan2() uses true pi, then
82 *
83 * asin(x) returns the exact asin(x) with error below about 2 ulps.
84 *
85 * In a test run with more than 1,024,000 random arguments on a VAX, the
86 * maximum observed error in ulps (units in the last place) was
87 * 1.99 ulps.
88 */
89
asin(x)90 double asin(x)
91 double x;
92 {
93 double s,t,copysign(),atan2(),sqrt(),one=1.0;
94 #if !defined(vax)&&!defined(tahoe)
95 if(x!=x) return(x); /* x is NaN */
96 #endif /* !defined(vax)&&!defined(tahoe) */
97 s=copysign(x,one);
98 if(s <= 0.5)
99 return(atan2(x,sqrt(one-x*x)));
100 else
101 { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
102
103 }
104
105 /* ACOS(X)
106 * RETURNS ARC COS OF X
107 * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
108 * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
109 *
110 * Required system supported functions:
111 * copysign(x,y)
112 * sqrt(x)
113 *
114 * Required kernel function:
115 * atan2(y,x)
116 *
117 * Method :
118 * ________
119 * / 1 - x
120 * acos(x) = 2*atan2( / -------- , 1 ) .
121 * \/ 1 + x
122 *
123 * Special cases:
124 * if x is NaN, return x itself;
125 * if |x|>1, return NaN.
126 *
127 * Accuracy:
128 * 1) If atan2() uses machine PI, then
129 *
130 * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
131 * and PI is the exact pi rounded to machine precision (see atan2 for
132 * details):
133 *
134 * in decimal:
135 * pi = 3.141592653589793 23846264338327 .....
136 * 53 bits PI = 3.141592653589793 115997963 ..... ,
137 * 56 bits PI = 3.141592653589793 227020265 ..... ,
138 *
139 * in hexadecimal:
140 * pi = 3.243F6A8885A308D313198A2E....
141 * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
142 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
143 *
144 * In a test run with more than 200,000 random arguments on a VAX, the
145 * maximum observed error in ulps (units in the last place) was
146 * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x)));
147 *
148 * 2) If atan2() uses true pi, then
149 *
150 * acos(x) returns the exact acos(x) with error below about 2 ulps.
151 *
152 * In a test run with more than 1,024,000 random arguments on a VAX, the
153 * maximum observed error in ulps (units in the last place) was
154 * 2.15 ulps.
155 */
156
acos(x)157 double acos(x)
158 double x;
159 {
160 double t,copysign(),atan2(),sqrt(),one=1.0;
161 #if !defined(vax)&&!defined(tahoe)
162 if(x!=x) return(x);
163 #endif /* !defined(vax)&&!defined(tahoe) */
164 if( x != -1.0)
165 t=atan2(sqrt((one-x)/(one+x)),one);
166 else
167 t=atan2(one,0.0); /* t = PI/2 */
168 return(t+t);
169 }
170