xref: /netbsd/lib/libm/src/s_sincos.c (revision 779858b8)
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 <sys/cdefs.h>
16 #if defined(LIBM_SCCS) && !defined(lint)
17 __RCSID("$NetBSD: s_sincos.c,v 1.5 2022/08/29 01:48:34 riastradh Exp $");
18 #endif
19 #if 0
20 __FBSDID("$FreeBSD: head/lib/msun/src/s_sincos.c 319047 2017-05-28 06:13:38Z mmel $");
21 #endif
22 
23 #include "namespace.h"
24 #include <float.h>
25 
26 #include "math.h"
27 // #define INLINE_REM_PIO2
28 #include "math_private.h"
29 // #include "e_rem_pio2.c"
30 #include "k_sincos.h"
31 
32 #ifdef __weak_alias
__weak_alias(sincos,_sincos)33 __weak_alias(sincos,_sincos)
34 #endif
35 
36 void
37 sincos(double x, double *sn, double *cs)
38 {
39 	double y[2];
40 	int32_t n, ix;
41 
42 	/* High word of x. */
43 	GET_HIGH_WORD(ix, x);
44 
45 	/* |x| ~< pi/4 */
46 	ix &= 0x7fffffff;
47 	if (ix <= 0x3fe921fb) {
48 		if (ix < 0x3e400000) {		/* |x| < 2**-27 */
49 			if ((int)x == 0) {	/* Generate inexact. */
50 				*sn = x;
51 				*cs = 1;
52 				return;
53 			}
54 		}
55 		__kernel_sincos(x, 0, 0, sn, cs);
56 		return;
57 	}
58 
59 	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
60 	if (ix >= 0x7ff00000) {
61 		*sn = x - x;
62 		*cs = x - x;
63 		return;
64 	}
65 
66 	/* Argument reduction. */
67 	n = __ieee754_rem_pio2(x, y);
68 
69 	switch(n & 3) {
70 	case 0:
71 		__kernel_sincos(y[0], y[1], 1, sn, cs);
72 		break;
73 	case 1:
74 		__kernel_sincos(y[0], y[1], 1, cs, sn);
75 		*cs = -*cs;
76 		break;
77 	case 2:
78 		__kernel_sincos(y[0], y[1], 1, sn, cs);
79 		*sn = -*sn;
80 		*cs = -*cs;
81 		break;
82 	default:
83 		__kernel_sincos(y[0], y[1], 1, cs, sn);
84 		*sn = -*sn;
85 	}
86 }
87 
88 #if !defined(__HAVE_LONG_DOUBLE)
89 __weak_alias(sincosl, sincos);
90 #endif
91