efficient sincos based on sin and cos

1 files changed, 66 insertions, 0 deletions

diff --git a/src/math/sincosl.c b/src/math/sincosl.c
new file mode 100644
index 00000000..378dc979
--- /dev/null
+++ b/src/math/sincosl.c
@@ -0,0 +1,66 @@
+#include "libm.h"
+
+#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
+void sincosl(long double x, long double *sin, long double *cos)
+{
+	double s, c;
+	sincos(x, &s, &c);
+	*sin = s;
+	*cos = c;
+}
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+#include "__rem_pio2l.h"
+
+void sincosl(long double x, long double *sin, long double *cos)
+{
+	union IEEEl2bits u;
+	int n;
+	long double y[2], s, c;
+
+	u.e = x;
+	u.bits.sign = 0;
+
+	/* x = +-0 or subnormal */
+	if (!u.bits.exp) {
+		*sin = x;
+		*cos = 1.0;
+		return;
+	}
+
+	/* x = nan or inf */
+	if (u.bits.exp == 0x7fff) {
+		*sin = *cos = x - x;
+		return;
+	}
+
+	/* |x| < pi/4 */
+	if (u.e < M_PI_4) {
+		*sin = __sinl(x, 0, 0);
+		*cos = __cosl(x, 0);
+		return;
+	}
+
+	n = __rem_pio2l(x, y);
+	s = __sinl(y[0], y[1], 1);
+	c = __cosl(y[0], y[1]);
+	switch (n & 3) {
+	case 0:
+		*sin = s;
+		*cos = c;
+		break;
+	case 1:
+		*sin = c;
+		*cos = -s;
+		break;
+	case 2:
+		*sin = -s;
+		*cos = -c;
+		break;
+	case 3:
+	default:
+		*sin = -c;
+		*cos = s;
+		break;
+	}
+}
+#endif