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