Implement sicos(3), sincosf(3) and sincosl(3). These functions are common

extensions and modern compilers (such as clang) will use them to optimize
separate calculations of sine and cosine.

ok tom@, patrick@, deraadt@, jmc@

7 files changed, 539 insertions, 0 deletions

diff --git a/lib/libm/src/k_sincos.h b/lib/libm/src/k_sincos.h
new file mode 100644
index 00000000000..796e72af89e
--- /dev/null
+++ b/lib/libm/src/k_sincos.h
@@ -0,0 +1,49 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ * 
+ * k_sin.c and k_cos.c merged by Steven G. Kargl.
+ */
+
+static const double
+S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
+S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
+S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
+S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
+S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
+S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
+
+static const double
+C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
+C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
+C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
+C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
+C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
+C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
+
+static inline void
+__kernel_sincos(double x, double y, int iy, double *sn, double *cs)
+{
+	double hz, r, v, w, z;
+
+	z = x * x;
+	w = z * z;
+	r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
+	v = z * x;
+
+	if (iy == 0)
+		*sn = x + v * (S1 + z * r);
+	else
+		*sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
+
+	r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
+	hz = z / 2;
+	w = 1 - hz;
+	*cs = w + (((1 - w) - hz) + (z * r - x * y));
+}
diff --git a/lib/libm/src/k_sincosf.h b/lib/libm/src/k_sincosf.h
new file mode 100644
index 00000000000..f031016b79d
--- /dev/null
+++ b/lib/libm/src/k_sincosf.h
@@ -0,0 +1,40 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * k_sinf.c and k_cosf.c merged by Steven G. Kargl.
+ */
+
+/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
+static const double
+S1 = -0x15555554cbac77.0p-55,	/* -0.166666666416265235595 */
+S2 =  0x111110896efbb2.0p-59,	/*  0.0083333293858894631756 */
+S3 = -0x1a00f9e2cae774.0p-65,	/* -0.000198393348360966317347 */
+S4 =  0x16cd878c3b46a7.0p-71;	/*  0.0000027183114939898219064 */
+
+/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
+static const double
+C0  = -0x1ffffffd0c5e81.0p-54,	/* -0.499999997251031003120 */
+C1  =  0x155553e1053a42.0p-57,	/*  0.0416666233237390631894 */
+C2  = -0x16c087e80f1e27.0p-62,	/* -0.00138867637746099294692 */
+C3  =  0x199342e0ee5069.0p-68;	/*  0.0000243904487962774090654 */
+
+static inline void
+__kernel_sincosdf(double x, float *sn, float *cs)
+{
+	double r, s, w, z;
+
+	z = x * x;
+	w = z * z;
+	r = S3 + z * S4;
+	s = z * x;
+	*sn = (x + s * (S1 + z * S2)) + s * w * r;
+	r = C2 + z * C3;
+	*cs = ((1 + z * C0) + w * C1) + (w * z) * r;
+}
diff --git a/lib/libm/src/ld128/k_sincosl.h b/lib/libm/src/ld128/k_sincosl.h
new file mode 100644
index 00000000000..55247bbd189
--- /dev/null
+++ b/lib/libm/src/ld128/k_sincosl.h
@@ -0,0 +1,68 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * k_sinl.c and k_cosl.c merged by Steven G. Kargl
+ */
+
+static const long double
+C1 =  0.04166666666666666666666666666666658424671L,
+C2 = -0.001388888888888888888888888888863490893732L,
+C3 =  0.00002480158730158730158730158600795304914210L,
+C4 = -0.2755731922398589065255474947078934284324e-6L,
+C5 =  0.2087675698786809897659225313136400793948e-8L,
+C6 = -0.1147074559772972315817149986812031204775e-10L,
+C7 =  0.4779477332386808976875457937252120293400e-13L,
+S1 = -0.16666666666666666666666666666666666606732416116558L,
+S2 =  0.0083333333333333333333333333333331135404851288270047L,
+S3 = -0.00019841269841269841269841269839935785325638310428717L,
+S4 =  0.27557319223985890652557316053039946268333231205686e-5L,
+S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
+S6 =  0.16059043836821614596571832194524392581082444805729e-9L,
+S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
+S8 =  0.28114572543451292625024967174638477283187397621303e-14L;
+
+static const double
+C8  = -0.1561920696721507929516718307820958119868e-15,
+C9  =  0.4110317413744594971475941557607804508039e-18,
+C10 = -0.8896592467191938803288521958313920156409e-21,
+C11 =  0.1601061435794535138244346256065192782581e-23,
+S9  = -0.82206352458348947812512122163446202498005154296863e-17,
+S10 =  0.19572940011906109418080609928334380560135358385256e-19,
+S11 = -0.38680813379701966970673724299207480965452616911420e-22,
+S12 =  0.64038150078671872796678569586315881020659912139412e-25;
+
+static inline void
+__kernel_sincosl(long double x, long double y, int iy, long double *sn, 
+    long double *cs)
+{
+	long double hz, r, v, w, z;
+
+	z = x * x;
+	v = z * x;
+	/*
+	 * XXX Replace Horner scheme with an algorithm suitable for CPUs
+	 * with more complex pipelines.
+	 */
+	r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
+	    z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
+
+	if (iy == 0)
+		*sn = x + v * (S1 + z * r);
+	else
+		*cs = x - ((z * (y / 2 - v * r) - y) - v * S1);
+
+	hz = z / 2;
+	w = 1 - hz;
+	r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + 
+	    z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
+
+	*cs =  w + (((1 - w) - hz) + (z * r - x * y));
+}
diff --git a/lib/libm/src/ld80/k_sincosl.h b/lib/libm/src/ld80/k_sincosl.h
new file mode 100644
index 00000000000..82421722f68
--- /dev/null
+++ b/lib/libm/src/ld80/k_sincosl.h
@@ -0,0 +1,69 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * k_sinl.c and k_cosl.c merged by Steven G. Kargl
+ */
+
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+C1hi = 0.041666666666666664,		/*  0x15555555555555.0p-57 */
+C1lo = 2.2598839032744733e-18,		/*  0x14d80000000000.0p-111 */
+S1hi = -0.16666666666666666,		/* -0x15555555555555.0p-55 */
+S1lo = -9.2563760475949941e-18;		/* -0x15580000000000.0p-109 */
+#define	S1	((long double)S1hi + S1lo)
+#define	C1	((long double)C1hi + C1lo)
+#else
+static const long double
+C1 =  0.0416666666666666666136L;	/*  0xaaaaaaaaaaaaaa9b.0p-68 */
+S1 = -0.166666666666666666671L,		/* -0xaaaaaaaaaaaaaaab.0p-66 */
+#endif
+
+static const double
+C2 = -0.0013888888888888874,		/* -0x16c16c16c16c10.0p-62 */
+C3 =  0.000024801587301571716,		/*  0x1a01a01a018e22.0p-68 */
+C4 = -0.00000027557319215507120,	/* -0x127e4fb7602f22.0p-74 */
+C5 =  0.0000000020876754400407278,	/*  0x11eed8caaeccf1.0p-81 */
+C6 = -1.1470297442401303e-11,		/* -0x19393412bd1529.0p-89 */
+C7 =  4.7383039476436467e-14,		/*  0x1aac9d9af5c43e.0p-97 */
+S2 =  0.0083333333333333332,		/*  0x11111111111111.0p-59 */
+S3 = -0.00019841269841269427,		/* -0x1a01a01a019f81.0p-65 */
+S4 =  0.0000027557319223597490,		/*  0x171de3a55560f7.0p-71 */
+S5 = -0.000000025052108218074604,	/* -0x1ae64564f16cad.0p-78 */
+S6 =  1.6059006598854211e-10,		/*  0x161242b90243b5.0p-85 */
+S7 = -7.6429779983024564e-13,		/* -0x1ae42ebd1b2e00.0p-93 */
+S8 =  2.6174587166648325e-15;		/*  0x179372ea0b3f64.0p-101 */
+
+static inline void
+__kernel_sincosl(long double x, long double y, int iy, long double *sn,
+    long double *cs)
+{
+	long double hz, r, v, w, z;
+
+	z = x * x;
+	v = z * x;
+	/*
+	 * XXX Replace Horner scheme with an algorithm suitable for CPUs
+	 * with more complex pipelines.
+	 */
+	r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * S8)))));
+
+	if (iy == 0)
+		*sn = x + v * (S1 + z * r);
+	else
+		*sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
+
+	hz = z / 2;
+	w = 1 - hz;
+	r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
+	    z * C7))))));
+	*cs = w + (((1 - w) - hz) + (z * r - x * y));
+}
diff --git a/lib/libm/src/s_sincos.c b/lib/libm/src/s_sincos.c
new file mode 100644
index 00000000000..754205c5626
--- /dev/null
+++ b/lib/libm/src/s_sincos.c
@@ -0,0 +1,72 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ * 
+ * s_sin.c and s_cos.c merged by Steven G. Kargl.  Descriptions of the
+ * algorithms are contained in the original files.
+ */
+
+#include <float.h>
+#include "math.h"
+
+#include "math_private.h"
+#include "k_sincos.h"
+
+void
+sincos(double x, double *sn, double *cs)
+{
+	double y[2];
+	int32_t n, ix;
+
+	/* High word of x. */
+	GET_HIGH_WORD(ix, x);
+
+	/* |x| ~< pi/4 */
+	ix &= 0x7fffffff;
+	if (ix <= 0x3fe921fb) {
+		if (ix < 0x3e400000) {		/* |x| < 2**-27 */
+			if ((int)x == 0) {	/* Generate inexact. */
+				*sn = x;
+				*cs = 1;
+				return;
+			}
+		}
+		__kernel_sincos(x, 0, 0, sn, cs);
+		return;
+	}
+
+	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
+	if (ix >= 0x7ff00000) {
+		*sn = x - x;
+		*cs = x - x;
+		return;
+	}
+
+	/* Argument reduction. */
+	n = __ieee754_rem_pio2(x, y);
+
+	switch(n & 3) {
+	case 0:
+		__kernel_sincos(y[0], y[1], 1, sn, cs);
+		break;
+	case 1:
+		__kernel_sincos(y[0], y[1], 1, cs, sn);
+		*cs = -*cs;
+		break;
+	case 2:
+		__kernel_sincos(y[0], y[1], 1, sn, cs);
+		*sn = -*sn;
+		*cs = -*cs;
+		break;
+	default:
+		__kernel_sincos(y[0], y[1], 1, cs, sn);
+		*sn = -*sn;
+	}
+}
+LDBL_MAYBE_CLONE(sincos);
diff --git a/lib/libm/src/s_sincosf.c b/lib/libm/src/s_sincosf.c
new file mode 100644
index 00000000000..a02aa92ef3b
--- /dev/null
+++ b/lib/libm/src/s_sincosf.c
@@ -0,0 +1,120 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* s_sincosf.c -- float version of s_sincos.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ * Merged s_sinf.c and s_cosf.c by Steven G. Kargl.
+ */
+
+#include <float.h>
+#include "math.h"
+
+#include "math_private.h"
+#include "k_sincosf.h"
+
+/* Small multiples of pi/2 rounded to double precision. */
+static const double
+p1pio2 = 1*M_PI_2,			/* 0x3FF921FB, 0x54442D18 */
+p2pio2 = 2*M_PI_2,			/* 0x400921FB, 0x54442D18 */
+p3pio2 = 3*M_PI_2,			/* 0x4012D97C, 0x7F3321D2 */
+p4pio2 = 4*M_PI_2;			/* 0x401921FB, 0x54442D18 */
+
+void
+sincosf(float x, float *sn, float *cs)
+{
+	float c, s;
+	float y[2];
+	int32_t n, hx, ix;
+
+	GET_FLOAT_WORD(hx, x);
+	ix = hx & 0x7fffffff;
+
+	if (ix <= 0x3f490fda) {		/* |x| ~<= pi/4 */
+		if (ix < 0x39800000) {	/* |x| < 2**-12 */
+			if ((int)x == 0) {
+				*sn = x;	/* x with inexact if x != 0 */
+				*cs = 1;
+				return;
+			}
+		}
+		__kernel_sincosdf(x, sn, cs);
+		return;
+	}
+
+	if (ix <= 0x407b53d1) {		/* |x| ~<= 5*pi/4 */
+		if (ix <= 0x4016cbe3) {	/* |x| ~<= 3pi/4 */
+			if (hx > 0) {
+				__kernel_sincosdf(x - p1pio2, cs, sn);
+				*cs = -*cs;
+			} else {
+				__kernel_sincosdf(x + p1pio2, cs, sn);
+				*sn = -*sn;
+			}
+		} else {
+			if (hx > 0)
+				__kernel_sincosdf(x - p2pio2, sn, cs);
+			else
+				__kernel_sincosdf(x + p2pio2, sn, cs);
+			*sn = -*sn;
+			*cs = -*cs;
+		}
+		return;
+	}
+
+	if (ix <= 0x40e231d5) {		/* |x| ~<= 9*pi/4 */
+		if (ix <= 0x40afeddf) {	/* |x| ~<= 7*pi/4 */
+			if (hx > 0) {
+				__kernel_sincosdf(x - p3pio2, cs, sn);
+				*sn = -*sn;
+			} else {
+				__kernel_sincosdf(x + p3pio2, cs, sn);
+				*cs = -*cs;
+			}
+		} else {
+			if (hx > 0)
+				__kernel_sincosdf(x - p4pio2, sn, cs);
+			else
+				__kernel_sincosdf(x + p4pio2, sn, cs);
+		}
+		return;
+	}
+
+	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
+	if (ix >= 0x7f800000) {
+		*sn = x - x;
+		*cs = x - x;
+		return;
+	}
+
+	/* Argument reduction. */
+	n = __ieee754_rem_pio2f(x, y);
+	s = __kernel_sinf(y[0], y[1], 1);
+	c = __kernel_cosf(y[0], y[1]);
+
+	switch(n & 3) {
+	case 0:
+		*sn = s;
+		*cs = c;
+		break;
+	case 1:
+		*sn = c;
+		*cs = -s;
+		break;
+	case 2:
+		*sn = -s;
+		*cs = -c;
+		break;
+	default:
+		*sn = -c;
+		*cs = s;
+	}
+}
diff --git a/lib/libm/src/s_sincosl.c b/lib/libm/src/s_sincosl.c
new file mode 100644
index 00000000000..e89d30c4b90
--- /dev/null
+++ b/lib/libm/src/s_sincosl.c
@@ -0,0 +1,121 @@
+/*-
+ * Copyright (c) 2007, 2010-2013 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * s_sinl.c and s_cosl.c merged by Steven G. Kargl.
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
+
+#if LDBL_MANT_DIG == 64
+#include "../ld80/k_sincosl.h"
+#define	NX	3
+#define	PREC	2
+#elif LDBL_MANT_DIG == 113
+#include "../ld128/k_sincosl.h"
+#define	NX	5
+#define	PREC	3
+#else
+#error "Unsupported long double format"
+#endif
+
+static const long double two24 = 1.67772160000000000000e+07L;
+
+void
+sincosl(long double x, long double *sn, long double *cs)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} z;
+	int i, e0;
+	double xd[NX], yd[PREC];
+	long double hi, lo;
+
+	z.e = x;
+	z.bits.ext_sign = 0;
+
+	/* Optimize the case where x is already within range. */
+	if (z.e < M_PI_4) {
+		/*
+		 * If x = +-0 or x is a subnormal number, then sin(x) = x and
+		 * cos(x) = 1.
+		 */
+		if (z.bits.ext_exp == 0) {
+			*sn = x;
+			*cs = 1;
+		} else
+			__kernel_sincosl(x, 0, 0, sn, cs);
+	}
+
+	/* If x = NaN or Inf, then sin(x) and cos(x) are NaN. */
+	if (z.bits.ext_exp == 32767) {
+		*sn = x - x;
+		*cs = x - x;
+	}
+
+	/* Split z.e into a 24-bit representation. */
+	e0 = ilogbl(z.e) - 23;
+	z.e = scalbnl(z.e, -e0);
+	for (i = 0; i < NX; i++) {
+		xd[i] = (double)((int32_t)z.e);
+		z.e = (z.e - xd[i]) * two24;
+	}
+
+	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
+	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
+
+#if PREC == 2
+	hi = (long double)yd[0] + yd[1];
+	lo = yd[1] - (hi - yd[0]);
+#else /* PREC == 3 */
+	long double t;
+	t = (long double)yd[2] + yd[1];
+	hi = t + yd[0];
+	lo = yd[0] - (hi - t);
+#endif
+
+	switch (e0 & 3) {
+	case 0:
+		__kernel_sincosl(hi, lo, 1, sn, cs);
+		break;
+	case 1:
+		__kernel_sincosl(hi, lo, 1, cs, sn);
+		*cs = -*cs;
+		break;
+	case 2:
+		__kernel_sincosl(hi, lo, 1, sn, cs);
+		*sn = -*sn;
+		*cs = -*cs;
+		break;
+	default:
+		__kernel_sincosl(hi, lo, 1, cs, sn);
+		*sn = -*sn;
+	}
+}