1 /*
2 * Copyright (c) 1985, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * %sccs.include.redist.c%
6 */
7
8 #ifndef lint
9 static char sccsid[] = "@(#)asincos.c 8.1 (Berkeley) 06/04/93";
10 #endif /* not lint */
11
12 /* ASIN(X)
13 * RETURNS ARC SINE OF X
14 * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
15 * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
16 *
17 * Required system supported functions:
18 * copysign(x,y)
19 * sqrt(x)
20 *
21 * Required kernel function:
22 * atan2(y,x)
23 *
24 * Method :
25 * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
26 * computed as follows
27 * 1-x*x if x < 0.5,
28 * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
29 *
30 * Special cases:
31 * if x is NaN, return x itself;
32 * if |x|>1, return NaN.
33 *
34 * Accuracy:
35 * 1) If atan2() uses machine PI, then
36 *
37 * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
38 * and PI is the exact pi rounded to machine precision (see atan2 for
39 * details):
40 *
41 * in decimal:
42 * pi = 3.141592653589793 23846264338327 .....
43 * 53 bits PI = 3.141592653589793 115997963 ..... ,
44 * 56 bits PI = 3.141592653589793 227020265 ..... ,
45 *
46 * in hexadecimal:
47 * pi = 3.243F6A8885A308D313198A2E....
48 * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
49 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
50 *
51 * In a test run with more than 200,000 random arguments on a VAX, the
52 * maximum observed error in ulps (units in the last place) was
53 * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x)));
54 *
55 * 2) If atan2() uses true pi, then
56 *
57 * asin(x) returns the exact asin(x) with error below about 2 ulps.
58 *
59 * In a test run with more than 1,024,000 random arguments on a VAX, the
60 * maximum observed error in ulps (units in the last place) was
61 * 1.99 ulps.
62 */
63
asin(x)64 double asin(x)
65 double x;
66 {
67 double s,t,copysign(),atan2(),sqrt(),one=1.0;
68 #if !defined(vax)&&!defined(tahoe)
69 if(x!=x) return(x); /* x is NaN */
70 #endif /* !defined(vax)&&!defined(tahoe) */
71 s=copysign(x,one);
72 if(s <= 0.5)
73 return(atan2(x,sqrt(one-x*x)));
74 else
75 { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
76
77 }
78
79 /* ACOS(X)
80 * RETURNS ARC COS OF X
81 * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
82 * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
83 *
84 * Required system supported functions:
85 * copysign(x,y)
86 * sqrt(x)
87 *
88 * Required kernel function:
89 * atan2(y,x)
90 *
91 * Method :
92 * ________
93 * / 1 - x
94 * acos(x) = 2*atan2( / -------- , 1 ) .
95 * \/ 1 + x
96 *
97 * Special cases:
98 * if x is NaN, return x itself;
99 * if |x|>1, return NaN.
100 *
101 * Accuracy:
102 * 1) If atan2() uses machine PI, then
103 *
104 * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
105 * and PI is the exact pi rounded to machine precision (see atan2 for
106 * details):
107 *
108 * in decimal:
109 * pi = 3.141592653589793 23846264338327 .....
110 * 53 bits PI = 3.141592653589793 115997963 ..... ,
111 * 56 bits PI = 3.141592653589793 227020265 ..... ,
112 *
113 * in hexadecimal:
114 * pi = 3.243F6A8885A308D313198A2E....
115 * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps
116 * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps
117 *
118 * In a test run with more than 200,000 random arguments on a VAX, the
119 * maximum observed error in ulps (units in the last place) was
120 * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x)));
121 *
122 * 2) If atan2() uses true pi, then
123 *
124 * acos(x) returns the exact acos(x) with error below about 2 ulps.
125 *
126 * In a test run with more than 1,024,000 random arguments on a VAX, the
127 * maximum observed error in ulps (units in the last place) was
128 * 2.15 ulps.
129 */
130
acos(x)131 double acos(x)
132 double x;
133 {
134 double t,copysign(),atan2(),sqrt(),one=1.0;
135 #if !defined(vax)&&!defined(tahoe)
136 if(x!=x) return(x);
137 #endif /* !defined(vax)&&!defined(tahoe) */
138 if( x != -1.0)
139 t=atan2(sqrt((one-x)/(one+x)),one);
140 else
141 t=atan2(one,0.0); /* t = PI/2 */
142 return(t+t);
143 }
144