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