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