xref: /freebsd/lib/msun/src/s_sincos.c (revision 8aac90f18aef7c9eea906c3ff9a001ca7b94f375)
1 /*-
2  * ====================================================
3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4  *
5  * Developed at SunPro, a Sun Microsystems, Inc. business.
6  * Permission to use, copy, modify, and distribute this
7  * software is freely granted, provided that this notice
8  * is preserved.
9  * ====================================================
10  *
11  * s_sin.c and s_cos.c merged by Steven G. Kargl.  Descriptions of the
12  * algorithms are contained in the original files.
13  */
14 
15 #include <float.h>
16 
17 #include "math.h"
18 #define INLINE_REM_PIO2
19 #include "math_private.h"
20 #include "e_rem_pio2.c"
21 #include "k_sincos.h"
22 
23 void
24 sincos(double x, double *sn, double *cs)
25 {
26 	double y[2];
27 	int32_t n, ix;
28 
29 	/* High word of x. */
30 	GET_HIGH_WORD(ix, x);
31 
32 	/* |x| ~< pi/4 */
33 	ix &= 0x7fffffff;
34 	if (ix <= 0x3fe921fb) {
35 		if (ix < 0x3e400000) {		/* |x| < 2**-27 */
36 			if ((int)x == 0) {	/* Generate inexact. */
37 				*sn = x;
38 				*cs = 1;
39 				return;
40 			}
41 		}
42 		__kernel_sincos(x, 0, 0, sn, cs);
43 		return;
44 	}
45 
46 	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
47 	if (ix >= 0x7ff00000) {
48 		*sn = x - x;
49 		*cs = x - x;
50 		return;
51 	}
52 
53 	/* Argument reduction. */
54 	n = __ieee754_rem_pio2(x, y);
55 
56 	switch(n & 3) {
57 	case 0:
58 		__kernel_sincos(y[0], y[1], 1, sn, cs);
59 		break;
60 	case 1:
61 		__kernel_sincos(y[0], y[1], 1, cs, sn);
62 		*cs = -*cs;
63 		break;
64 	case 2:
65 		__kernel_sincos(y[0], y[1], 1, sn, cs);
66 		*sn = -*sn;
67 		*cs = -*cs;
68 		break;
69 	default:
70 		__kernel_sincos(y[0], y[1], 1, cs, sn);
71 		*sn = -*sn;
72 	}
73 }
74 
75 #if (LDBL_MANT_DIG == 53)
76 __weak_reference(sincos, sincosl);
77 #endif
78