math.h: make __FLOAT_BITS and __DOUBLE_BITS C89
[musl] / src / math / sincosl.c
1 #define _GNU_SOURCE
2 #include "libm.h"
3
4 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
5 void sincosl(long double x, long double *sin, long double *cos)
6 {
7         sincos(x, (double *)sin, (double *)cos);
8 }
9 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
10 void sincosl(long double x, long double *sin, long double *cos)
11 {
12         union ldshape u = {x};
13         unsigned n;
14         long double y[2], s, c;
15
16         u.i.se &= 0x7fff;
17         if (u.i.se == 0x7fff) {
18                 *sin = *cos = x - x;
19                 return;
20         }
21         if (u.f < M_PI_4) {
22                 if (u.i.se < 0x3fff - LDBL_MANT_DIG) {
23                         /* raise underflow if subnormal */
24                         if (u.i.se == 0) FORCE_EVAL(x*0x1p-120f);
25                         *sin = x;
26                         /* raise inexact if x!=0 */
27                         *cos = 1.0 + x;
28                         return;
29                 }
30                 *sin = __sinl(x, 0, 0);
31                 *cos = __cosl(x, 0);
32                 return;
33         }
34         n = __rem_pio2l(x, y);
35         s = __sinl(y[0], y[1], 1);
36         c = __cosl(y[0], y[1]);
37         switch (n & 3) {
38         case 0:
39                 *sin = s;
40                 *cos = c;
41                 break;
42         case 1:
43                 *sin = c;
44                 *cos = -s;
45                 break;
46         case 2:
47                 *sin = -s;
48                 *cos = -c;
49                 break;
50         case 3:
51         default:
52                 *sin = -c;
53                 *cos = s;
54                 break;
55         }
56 }
57 #endif