xref: /csrg-svn/lib/libm/common/sincos.c (revision 31931)
1 /*
2  * Copyright (c) 1987 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[] = "@(#)sincos.c	1.1	1.1 (ucb.elefunt) 07/24/87";
16 #endif	/* not lint */
17 
18 #include "trig.h"
19 double
20 sin(x)
21 double x;
22 {
23 	double a,c,z;
24 
25         if(!finite(x))		/* sin(NaN) and sin(INF) must be NaN */
26 		return x-x;
27 	x=drem(x,PI2);		/* reduce x into [-PI,PI] */
28 	a=copysign(x,one);
29 	if (a >= PIo4) {
30 		if(a >= PI3o4)		/* ... in [3PI/4,PI] */
31 			x = copysign((a = PI-a),x);
32 		else {			/* ... in [PI/4,3PI/4]  */
33 			a = PIo2-a;		/* rtn. sign(x)*C(PI/2-|x|) */
34 			z = a*a;
35 			c = cos__C(z);
36 			z *= half;
37 			a = (z >= thresh ? half-((z-half)-c) : one-(z-c));
38 			return copysign(a,x);
39 		}
40 	}
41 
42 	if (a < small) {		/* rtn. S(x) */
43 		big+a;
44 		return x;
45 	}
46 	return x+x*sin__S(x*x);
47 }
48 
49 double
50 cos(x)
51 double x;
52 {
53 	double a,c,z,s = 1.0;
54 
55 	if(!finite(x))		/* cos(NaN) and cos(INF) must be NaN */
56 		return x-x;
57 	x=drem(x,PI2);		/* reduce x into [-PI,PI] */
58 	a=copysign(x,one);
59 	if (a >= PIo4) {
60 		if (a >= PI3o4) {	/* ... in [3PI/4,PI] */
61 			a = PI-a;
62 			s = negone;
63 		}
64 		else {			/* ... in [PI/4,3PI/4] */
65 			a = PIo2-a;
66 			return a+a*sin__S(a*a);	/* rtn. S(PI/2-|x|) */
67 		}
68 	}
69 	if (a < small) {
70 		big+a;
71 		return s;		/* rtn. s*C(a) */
72 	}
73 	z = a*a;
74 	c = cos__C(z);
75 	z *= half;
76 	a = (z >= thresh ? half-((z-half)-c) : one-(z-c));
77 	return copysign(a,s);
78 }
79