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