The diff below adds sincos(3) to our libm. This function isn't
standardized but is found on Linux, MacOS X and FreeBSD. It is
essentially an optimization. Calculating both sin and cos of the same
numbers is a fairly common operation. By calculating both in the same
function, argument reduction has to be done only once. And argument
reduction is a fairly expensive operation.
Modern compilers know about this function and will replace separate
sin(3) and cos(3) calls with a sincos(3) call if they know the latter
is available. In principle our clang knows about this as well, but it
"knows" this function isn't available on OpenBSD. Except on armv7,
because we claim to be armv7-unknown-openbsd6.3-gnueabi and clang
thinks that gnueabi implies a "GNU" environment that has sincos(3).
Teaching it that OpenBSD is "special" isn't hard, but fixing this
properly opens a can of worms. Adding sincos(3) to libm is easier and
happens to be the route chosen by FreeBSD.
Diff below is adapted from FreeBSD. Still needs a bit of testing and
a man page.
Thoughts?
Index: include/math.h
===================================================================
RCS file: /cvs/src/include/math.h,v
retrieving revision 1.35
diff -u -p -r1.35 math.h
--- include/math.h 17 Mar 2016 20:55:35 -0000 1.35
+++ include/math.h 2 Mar 2018 21:00:20 -0000
@@ -273,6 +273,8 @@ int finite(double);
double gamma_r(double, int *);
double lgamma_r(double, int *);
+void sincos(double, double *, double *);
+
/*
* IEEE Test Vector
*/
@@ -382,6 +384,8 @@ int isnanf(float);
float gammaf_r(float, int *);
float lgammaf_r(float, int *);
+void sincosf(float, float *, float *);
+
/*
* Float version of IEEE Test Vector
*/
@@ -459,6 +463,13 @@ long double fminl(long double, long doub
long double fmal(long double, long double, long double);
#endif /* __ISO_C_VISIBLE >= 1999 */
+
+/*
+ * Long double versions of BSD math library entry points
+ */
+#if __BSD_VISIBLE
+void sincosl(long double, long double *, long double *);
+#endif
/*
* Library implementation
Index: lib/libm/Makefile
===================================================================
RCS file: /cvs/src/lib/libm/Makefile,v
retrieving revision 1.115
diff -u -p -r1.115 Makefile
--- lib/libm/Makefile 11 Jan 2017 13:38:13 -0000 1.115
+++ lib/libm/Makefile 2 Mar 2018 21:00:20 -0000
@@ -101,8 +101,8 @@ COMMON_SRCS = b_exp__D.c b_log__D.c b_tg
s_nextafterf.c s_nexttowardf.c s_remquo.c s_remquof.c s_rint.c \
s_rintf.c \
s_scalbn.c s_scalbnf.c s_significand.c \
- s_significandf.c \
- s_sin.c s_sinf.c s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_trunc.c \
+ s_significandf.c s_sin.c s_sinf.c s_sincos.c s_sincosf.c \
+ s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_trunc.c \
s_truncf.c w_drem.c w_dremf.c w_gamma.c w_gamma_r.c w_gammaf.c \
w_gammaf_r.c w_lgamma.c w_lgammaf.c
@@ -119,7 +119,7 @@ LONG_SRCS = e_acoshl.c e_acosl.c e_asinl
s_floorl.c s_fmal.c s_fmaxl.c s_fminl.c s_frexpl.c s_ilogbl.c \
s_llrintl.c s_llroundl.c s_log1pl.c s_logbl.c s_lrintl.c \
s_lroundl.c s_modfl.c s_nanl.c s_nextafterl.c s_nexttoward.c \
- s_remquol.c s_rintl.c s_roundl.c s_scalbnl.c s_sinl.c \
+ s_remquol.c s_rintl.c s_roundl.c s_scalbnl.c s_sinl.c s_sincosl.c \
s_tanhl.c s_tanl.c s_truncl.c
# math routines that are completely MI
Index: lib/libm/Symbols.map
===================================================================
RCS file: /cvs/src/lib/libm/Symbols.map,v
retrieving revision 1.1
diff -u -p -r1.1 Symbols.map
--- lib/libm/Symbols.map 12 Sep 2016 19:47:01 -0000 1.1
+++ lib/libm/Symbols.map 2 Mar 2018 21:00:20 -0000
@@ -255,6 +255,9 @@
significand;
significandf;
sin;
+ sincos;
+ sincosf;
+ sincosl;
sinf;
sinh;
sinhf;
Index: lib/libm/shlib_version
===================================================================
RCS file: /cvs/src/lib/libm/shlib_version,v
retrieving revision 1.20
diff -u -p -r1.20 shlib_version
--- lib/libm/shlib_version 12 Sep 2016 19:47:01 -0000 1.20
+++ lib/libm/shlib_version 2 Mar 2018 21:00:20 -0000
@@ -1,2 +1,2 @@
major=10
-minor=0
+minor=1
Index: lib/libm/hidden/math.h
===================================================================
RCS file: /cvs/src/lib/libm/hidden/math.h,v
retrieving revision 1.1
diff -u -p -r1.1 math.h
--- lib/libm/hidden/math.h 12 Sep 2016 19:47:02 -0000 1.1
+++ lib/libm/hidden/math.h 2 Mar 2018 21:00:20 -0000
@@ -172,6 +172,9 @@ PROTO_NORMAL(scalbnl);
PROTO_DEPRECATED(significand);
PROTO_DEPRECATED(significandf);
PROTO_NORMAL(sin);
+PROTO_DEPRECATED(sincos);
+PROTO_DEPRECATED(sincosf);
+PROTO_DEPRECATED(sincosl);
PROTO_NORMAL(sinf);
PROTO_NORMAL(sinh);
PROTO_NORMAL(sinhf);
Index: lib/libm/src/k_sincos.h
===================================================================
RCS file: lib/libm/src/k_sincos.h
diff -N lib/libm/src/k_sincos.h
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/k_sincos.h 2 Mar 2018 21:00:20 -0000
@@ -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));
+}
Index: lib/libm/src/k_sincosf.h
===================================================================
RCS file: lib/libm/src/k_sincosf.h
diff -N lib/libm/src/k_sincosf.h
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/k_sincosf.h 2 Mar 2018 21:00:20 -0000
@@ -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;
+}
Index: lib/libm/src/s_sincos.c
===================================================================
RCS file: lib/libm/src/s_sincos.c
diff -N lib/libm/src/s_sincos.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/s_sincos.c 2 Mar 2018 21:00:20 -0000
@@ -0,0 +1,73 @@
+/*-
+ * ====================================================
+ * 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;
+ }
+}
+DEF_NONSTD(sincos);
+LDBL_MAYBE_CLONE(sincos);
Index: lib/libm/src/s_sincosf.c
===================================================================
RCS file: lib/libm/src/s_sincosf.c
diff -N lib/libm/src/s_sincosf.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/s_sincosf.c 2 Mar 2018 21:00:20 -0000
@@ -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, i...@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;
+ double y;
+ 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);
+ __kernel_sincosdf(y, &s, &c);
+
+ 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;
+ }
+}
+DEF_NONSTD(sincosf);
Index: lib/libm/src/s_sincosl.c
===================================================================
RCS file: lib/libm/src/s_sincosl.c
diff -N lib/libm/src/s_sincosl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/s_sincosl.c 2 Mar 2018 21:00:20 -0000
@@ -0,0 +1,122 @@
+/*-
+ * 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;
+ }
+}
+DEF_NONSTD(sincosl);
Index: lib/libm/src/ld128/k_sincosl.h
===================================================================
RCS file: lib/libm/src/ld128/k_sincosl.h
diff -N lib/libm/src/ld128/k_sincosl.h
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/ld128/k_sincosl.h 2 Mar 2018 21:00:20 -0000
@@ -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));
+}
Index: lib/libm/src/ld80/k_sincosl.h
===================================================================
RCS file: lib/libm/src/ld80/k_sincosl.h
diff -N lib/libm/src/ld80/k_sincosl.h
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/ld80/k_sincosl.h 2 Mar 2018 21:00:20 -0000
@@ -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));
+}
