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) 08/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
asin(x)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
acos(x)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 }
150