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