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