math: fix __fpclassifyl(-0.0) for IEEE binary128
[oweals/musl.git] / src / math / sincos.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunPro, a Sun Microsystems, Inc. business.
7  * Permission to use, copy, modify, and distribute this
8  * software is freely granted, provided that this notice
9  * is preserved.
10  * ====================================================
11  */
12
13 #define _GNU_SOURCE
14 #include "libm.h"
15
16 void sincos(double x, double *sin, double *cos)
17 {
18         double y[2], s, c;
19         uint32_t ix;
20         unsigned n;
21
22         GET_HIGH_WORD(ix, x);
23         ix &= 0x7fffffff;
24
25         /* |x| ~< pi/4 */
26         if (ix <= 0x3fe921fb) {
27                 /* if |x| < 2**-27 * sqrt(2) */
28                 if (ix < 0x3e46a09e) {
29                         /* raise inexact if x!=0 and underflow if subnormal */
30                         FORCE_EVAL(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f);
31                         *sin = x;
32                         *cos = 1.0;
33                         return;
34                 }
35                 *sin = __sin(x, 0.0, 0);
36                 *cos = __cos(x, 0.0);
37                 return;
38         }
39
40         /* sincos(Inf or NaN) is NaN */
41         if (ix >= 0x7ff00000) {
42                 *sin = *cos = x - x;
43                 return;
44         }
45
46         /* argument reduction needed */
47         n = __rem_pio2(x, y);
48         s = __sin(y[0], y[1], 1);
49         c = __cos(y[0], y[1]);
50         switch (n&3) {
51         case 0:
52                 *sin = s;
53                 *cos = c;
54                 break;
55         case 1:
56                 *sin = c;
57                 *cos = -s;
58                 break;
59         case 2:
60                 *sin = -s;
61                 *cos = -c;
62                 break;
63         case 3:
64         default:
65                 *sin = -c;
66                 *cos = s;
67                 break;
68         }
69 }