xref: /dragonfly/contrib/openbsd_libm/src/s_sincos.c (revision 1b963492d9a42a3a322a1de102e0cb457dd626f1)
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 #include "math.h"
17 
18 #include "math_private.h"
19 #include "k_sincos.h"
20 
21 void
sincos(double x,double * sn,double * cs)22 sincos(double x, double *sn, double *cs)
23 {
24           double y[2];
25           int32_t n, ix;
26 
27           /* High word of x. */
28           GET_HIGH_WORD(ix, x);
29 
30           /* |x| ~< pi/4 */
31           ix &= 0x7fffffff;
32           if (ix <= 0x3fe921fb) {
33                     if (ix < 0x3e400000) {                  /* |x| < 2**-27 */
34                               if ((int)x == 0) {  /* Generate inexact. */
35                                         *sn = x;
36                                         *cs = 1;
37                                         return;
38                               }
39                     }
40                     __kernel_sincos(x, 0, 0, sn, cs);
41                     return;
42           }
43 
44           /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
45           if (ix >= 0x7ff00000) {
46                     *sn = x - x;
47                     *cs = x - x;
48                     return;
49           }
50 
51           /* Argument reduction. */
52           n = __ieee754_rem_pio2(x, y);
53 
54           switch(n & 3) {
55           case 0:
56                     __kernel_sincos(y[0], y[1], 1, sn, cs);
57                     break;
58           case 1:
59                     __kernel_sincos(y[0], y[1], 1, cs, sn);
60                     *cs = -*cs;
61                     break;
62           case 2:
63                     __kernel_sincos(y[0], y[1], 1, sn, cs);
64                     *sn = -*sn;
65                     *cs = -*cs;
66                     break;
67           default:
68                     __kernel_sincos(y[0], y[1], 1, cs, sn);
69                     *sn = -*sn;
70           }
71 }
72