1 /* @(#)s_cos.c 5.1 93/09/24 */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunPro, a Sun Microsystems, Inc. business.
7  * Permission to use, copy, modify, and distribute this
8  * software is freely granted, provided that this notice
9  * is preserved.
10  * ====================================================
11  */
12 
13 #include <sys/cdefs.h>
14 __FBSDID("$FreeBSD$");
15 
16 /* cos(x)
17  * Return cosine function of x.
18  *
19  * kernel function:
20  *	__kernel_sin		... sine function on [-pi/4,pi/4]
21  *	__kernel_cos		... cosine function on [-pi/4,pi/4]
22  *	__ieee754_rem_pio2	... argument reduction routine
23  *
24  * Method.
25  *      Let S,C and T denote the sin, cos and tan respectively on
26  *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
27  *	in [-pi/4 , +pi/4], and let n = k mod 4.
28  *	We have
29  *
30  *          n        sin(x)      cos(x)        tan(x)
31  *     ----------------------------------------------------------
32  *	    0	       S	   C		 T
33  *	    1	       C	  -S		-1/T
34  *	    2	      -S	  -C		 T
35  *	    3	      -C	   S		-1/T
36  *     ----------------------------------------------------------
37  *
38  * Special cases:
39  *      Let trig be any of sin, cos, or tan.
40  *      trig(+-INF)  is NaN, with signals;
41  *      trig(NaN)    is that NaN;
42  *
43  * Accuracy:
44  *	TRIG(x) returns trig(x) nearly rounded
45  */
46 
47 #include <float.h>
48 
49 #include "math.h"
50 #define INLINE_REM_PIO2
51 #include "math_private.h"
52 #include "e_rem_pio2.c"
53 
54 double
cos(double x)55 cos(double x)
56 {
57 	double y[2],z=0.0;
58 	int32_t n, ix;
59 
60     /* High word of x. */
61 	GET_HIGH_WORD(ix,x);
62 
63     /* |x| ~< pi/4 */
64 	ix &= 0x7fffffff;
65 	if(ix <= 0x3fe921fb) {
66 	    if(ix<0x3e46a09e)			/* if x < 2**-27 * sqrt(2) */
67 		if(((int)x)==0) return 1.0;	/* generate inexact */
68 	    return __kernel_cos(x,z);
69 	}
70 
71     /* cos(Inf or NaN) is NaN */
72 	else if (ix>=0x7ff00000) return x-x;
73 
74     /* argument reduction needed */
75 	else {
76 	    n = __ieee754_rem_pio2(x,y);
77 	    switch(n&3) {
78 		case 0: return  __kernel_cos(y[0],y[1]);
79 		case 1: return -__kernel_sin(y[0],y[1],1);
80 		case 2: return -__kernel_cos(y[0],y[1]);
81 		default:
82 		        return  __kernel_sin(y[0],y[1],1);
83 	    }
84 	}
85 }
86 
87 #if (LDBL_MANT_DIG == 53)
88 __weak_reference(cos, cosl);
89 #endif
90