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