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 
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 
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