This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: [PATCH] x86-64: Add sincosf with vector FMA
- From: "H.J. Lu" <hjl dot tools at gmail dot com>
- To: GNU C Library <libc-alpha at sourceware dot org>
- Date: Mon, 8 Jan 2018 08:03:10 -0800
- Subject: Re: [PATCH] x86-64: Add sincosf with vector FMA
- Authentication-results: sourceware.org; auth=none
- References: <20171220220054.GA16094@intel.com> <20171221003115.GA2338@intel.com> <CAMe9rOqPzZsohDA3jGv_hY7rS-QuWSkxLXYpDGxV5pP5nxcenQ@mail.gmail.com>
On Mon, Jan 8, 2018 at 4:43 AM, H.J. Lu <hjl.tools@gmail.com> wrote:
> On Wed, Dec 20, 2017 at 4:31 PM, H.J. Lu <hongjiu.lu@intel.com> wrote:
>> On Wed, Dec 20, 2017 at 02:00:54PM -0800, H.J. Lu wrote:
>>> Since the x86-64 assembly version of sincosf is higly optimized with
>>> vector instructions, there isn't much room for improvement. However
>>> s_sincosf.c written in C with vector math and intrinsics can be
>>> optimized by GCC with FMA.
>>>
>>> On Skylake, bench-sincosf reports performance improvement:
>>>
>>> Assembly FMA improvement
>>> max 104.042 106.614 -2%
>>> min 9.426 8.586 10%
>>> mean 20.6209 18.803 10%
>>>
>>> Any coments?
>>>
>>> H.J.
>>> * sysdeps/x86_64/fpu/multiarch/Makefile (libm-sysdep_routines):
>>> Add s_sincosf-sse2 and s_sincosf-fma.
>>> (CFLAGS-s_sincosf-fma.c): New.
>>> * sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: New file.
>>> * sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S: Likewise.
>>> * sysdeps/x86_64/fpu/multiarch/s_sincosf.c: Likewise.
>>> * sysdeps/x86_64/fpu/s_sincosf.S: Don't add alias if
>>> __sincosf is defined.
>>
>> Updated patch without typedef of __v2df and __v4sf which have been
>> provided in <x86intrin.h>. Tested with GCC 4.9/5/6/7 on x86-64.
>>
>> H.J.
>> ----
>> Since the x86-64 assembly version of sincosf is higly optimized with
>> vector instructions, there isn't much room for improvement. However
>> s_sincosf.c written in C with vector math and intrinsics can be
>> optimized by GCC with FMA.
>>
>> On Skylake, bench-sincosf reports performance improvement:
>>
>> Assembly FMA improvement
>> max 104.042 106.614 -2%
>> min 9.426 8.586 10%
>> mean 20.6209 18.803 10%
>>
>> * sysdeps/x86_64/fpu/multiarch/Makefile (libm-sysdep_routines):
>> Add s_sincosf-sse2 and s_sincosf-fma.
>> (CFLAGS-s_sincosf-fma.c): New.
>> * sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: New file.
>> * sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S: Likewise.
>> * sysdeps/x86_64/fpu/multiarch/s_sincosf.c: Likewise.
>> * sysdeps/x86_64/fpu/s_sincosf.S: Don't add alias if
>> __sincosf is defined.
>
> If there are no objections, I am checking it in.
This is the patch I am checking. I vectoried 2 double -> float conversions.
On Skylake, now I got bench-sincosf reports performance improvement:
Assembly FMA improvement
max 104.042 101.008 3%
min 9.426 8.586 10%
mean 20.6209 18.2238 13%
--
H.J.
From 21b5bca0c20b55525adb1fe8373a61e87e441b5f Mon Sep 17 00:00:00 2001
From: "H.J. Lu" <hjl.tools@gmail.com>
Date: Wed, 20 Dec 2017 11:06:35 -0800
Subject: [PATCH] x86-64: Add sincosf with vector FMA
Since the x86-64 assembly version of sincosf is higly optimized with
vector instructions, there isn't much room for improvement. However
s_sincosf.c written in C with vector math and intrinsics can be
optimized by GCC with FMA.
On Skylake, bench-sincosf reports performance improvement:
Assembly FMA improvement
max 104.042 101.008 3%
min 9.426 8.586 10%
mean 20.6209 18.2238 13%
* sysdeps/x86_64/fpu/multiarch/Makefile (libm-sysdep_routines):
Add s_sincosf-sse2 and s_sincosf-fma.
(CFLAGS-s_sincosf-fma.c): New.
* sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: New file.
* sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S: Likewise.
* sysdeps/x86_64/fpu/multiarch/s_sincosf.c: Likewise.
* sysdeps/x86_64/fpu/s_sincosf.S: Don't add alias if
__sincosf is defined.
---
NEWS | 4 +-
sysdeps/x86_64/fpu/multiarch/Makefile | 5 +-
sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c | 240 ++++++++++++++++++++++++++
sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S | 2 +
sysdeps/x86_64/fpu/multiarch/s_sincosf.c | 28 +++
sysdeps/x86_64/fpu/s_sincosf.S | 2 +
6 files changed, 277 insertions(+), 4 deletions(-)
create mode 100644 sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c
create mode 100644 sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S
create mode 100644 sysdeps/x86_64/fpu/multiarch/s_sincosf.c
diff --git a/NEWS b/NEWS
index 0d6f826dcc..75bf46776a 100644
--- a/NEWS
+++ b/NEWS
@@ -19,8 +19,8 @@ Major new features:
supported on i386, x86_64, x32 and aarch64.
* Optimized x86-64 asin, atan2, exp, expf, log, pow, atan, sin, cosf,
- sinf and tan with FMA, contributed by Arjan van de Ven and H.J. Lu
- from Intel.
+ sinf, sincosf and tan with FMA, contributed by Arjan van de Ven and
+ H.J. Lu from Intel.
* Optimized x86-64 trunc and truncf for processors with SSE4.1.
diff --git a/sysdeps/x86_64/fpu/multiarch/Makefile b/sysdeps/x86_64/fpu/multiarch/Makefile
index 0825340c0c..9a89bfc286 100644
--- a/sysdeps/x86_64/fpu/multiarch/Makefile
+++ b/sysdeps/x86_64/fpu/multiarch/Makefile
@@ -37,10 +37,10 @@ CFLAGS-slowpow-fma.c = -mfma -mavx2
CFLAGS-s_sin-fma.c = -mfma -mavx2
CFLAGS-s_tan-fma.c = -mfma -mavx2
-libm-sysdep_routines += s_sinf-sse2 s_cosf-sse2
+libm-sysdep_routines += s_sinf-sse2 s_cosf-sse2 s_sincosf-sse2
libm-sysdep_routines += e_exp2f-fma e_expf-fma e_log2f-fma e_logf-fma \
- e_powf-fma s_sinf-fma s_cosf-fma
+ e_powf-fma s_sinf-fma s_cosf-fma s_sincosf-fma
CFLAGS-e_exp2f-fma.c = -mfma -mavx2
CFLAGS-e_expf-fma.c = -mfma -mavx2
@@ -49,6 +49,7 @@ CFLAGS-e_logf-fma.c = -mfma -mavx2
CFLAGS-e_powf-fma.c = -mfma -mavx2
CFLAGS-s_sinf-fma.c = -mfma -mavx2
CFLAGS-s_cosf-fma.c = -mfma -mavx2
+CFLAGS-s_sincosf-fma.c = -mfma -mavx2
libm-sysdep_routines += e_exp-fma4 e_log-fma4 e_pow-fma4 s_atan-fma4 \
e_asin-fma4 e_atan2-fma4 s_sin-fma4 s_tan-fma4 \
diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c
new file mode 100644
index 0000000000..64abe7abca
--- /dev/null
+++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c
@@ -0,0 +1,240 @@
+/* Compute sine and cosine of argument optimized with vector.
+ Copyright (C) 2017 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <errno.h>
+#include <math.h>
+#include <math_private.h>
+#include <x86intrin.h>
+#include <libm-alias-float.h>
+#include "s_sincosf.h"
+
+#define SINCOSF __sincosf_fma
+
+#ifndef SINCOSF
+# define SINCOSF_FUNC __sincosf
+#else
+# define SINCOSF_FUNC SINCOSF
+#endif
+
+/* Chebyshev constants for sin and cos, range -PI/4 - PI/4. */
+static const __v2df V0 = { -0x1.5555555551cd9p-3, -0x1.ffffffffe98aep-2};
+static const __v2df V1 = { 0x1.1111110c2688bp-7, 0x1.55555545c50c7p-5 };
+static const __v2df V2 = { -0x1.a019f8b4bd1f9p-13, -0x1.6c16b348b6874p-10 };
+static const __v2df V3 = { 0x1.71d7264e6b5b4p-19, 0x1.a00eb9ac43ccp-16 };
+static const __v2df V4 = { -0x1.a947e1674b58ap-26, -0x1.23c97dd8844d7p-22 };
+
+/* Chebyshev constants for sin and cos, range 2^-27 - 2^-5. */
+static const __v2df VC0 = { -0x1.555555543d49dp-3, -0x1.fffffff5cc6fdp-2 };
+static const __v2df VC1 = { 0x1.110f475cec8c5p-7, 0x1.55514b178dac5p-5 };
+
+static const __v2df v2ones = { 1.0, 1.0 };
+
+/* Compute the sine and cosine values using Chebyshev polynomials where
+ THETA is the range reduced absolute value of the input
+ and it is less than Pi/4,
+ N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide
+ whether a sine or cosine approximation is more accurate and
+ SIGNBIT is used to add the correct sign after the Chebyshev
+ polynomial is computed. */
+static void
+reduced_sincos (const double theta, const unsigned int n,
+ const unsigned int signbit, float *sinx, float *cosx)
+{
+ __v2df v2x, v2sx, v2cx;
+ const __v2df v2theta = { theta, theta };
+ const __v2df v2theta2 = v2theta * v2theta;
+ /* Here sinf() and cosf() are calculated using sin Chebyshev polynomial:
+ x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */
+ v2x = V3 + v2theta2 * V4; /* S3+x^2*S4. */
+ v2x = V2 + v2theta2 * v2x; /* S2+x^2*(S3+x^2*S4). */
+ v2x = V1 + v2theta2 * v2x; /* S1+x^2*(S2+x^2*(S3+x^2*S4)). */
+ v2x = V0 + v2theta2 * v2x; /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))). */
+ v2x = v2theta2 * v2x;
+ v2cx = v2ones + v2x;
+ v2sx = v2theta + v2theta * v2x;
+ /* We are operating on |x|, so we need to add back the original
+ signbit for sinf. */
+ /* Determine positive or negative primary interval. */
+ /* Are we in the primary interval of sin or cos? */
+ if ((n & 2) == 0)
+ {
+ const __v2df v2sign =
+ {
+ ones[((n >> 2) & 1) ^ signbit],
+ ones[((n + 2) >> 2) & 1]
+ };
+ v2cx[0] = v2sx[0];
+ v2cx *= v2sign;
+ __v4sf v4sx = _mm_cvtpd_ps (v2cx);
+ *sinx = v4sx[0];
+ *cosx = v4sx[1];
+ }
+ else
+ {
+ const __v2df v2sign =
+ {
+ ones[((n + 2) >> 2) & 1],
+ ones[((n >> 2) & 1) ^ signbit]
+ };
+ v2cx[0] = v2sx[0];
+ v2cx *= v2sign;
+ __v4sf v4sx = _mm_cvtpd_ps (v2cx);
+ *sinx = v4sx[1];
+ *cosx = v4sx[0];
+ }
+}
+
+void
+SINCOSF_FUNC (float x, float *sinx, float *cosx)
+{
+ double theta = x;
+ double abstheta = fabs (theta);
+ uint32_t ix, xi;
+ GET_FLOAT_WORD (xi, x);
+ /* |x| */
+ ix = xi & 0x7fffffff;
+ /* If |x|< Pi/4. */
+ if (ix < 0x3f490fdb)
+ {
+ if (ix >= 0x3d000000) /* |x| >= 2^-5. */
+ {
+ __v2df v2x, v2sx, v2cx;
+ const __v2df v2theta = { theta, theta };
+ const __v2df v2theta2 = v2theta * v2theta;
+ /* Chebyshev polynomial of the form for sin and cos. */
+ v2x = V3 + v2theta2 * V4;
+ v2x = V2 + v2theta2 * v2x;
+ v2x = V1 + v2theta2 * v2x;
+ v2x = V0 + v2theta2 * v2x;
+ v2x = v2theta2 * v2x;
+ v2cx = v2ones + v2x;
+ v2sx = v2theta + v2theta * v2x;
+ v2cx[0] = v2sx[0];
+ __v4sf v4sx = _mm_cvtpd_ps (v2cx);
+ *sinx = v4sx[0];
+ *cosx = v4sx[1];
+ }
+ else if (ix >= 0x32000000) /* |x| >= 2^-27. */
+ {
+ /* A simpler Chebyshev approximation is close enough for this range:
+ for sin: x+x^3*(SS0+x^2*SS1)
+ for cos: 1.0+x^2*(CC0+x^3*CC1). */
+ __v2df v2x, v2sx, v2cx;
+ const __v2df v2theta = { theta, theta };
+ const __v2df v2theta2 = v2theta * v2theta;
+ v2x = VC0 + v2theta * v2theta2 * VC1;
+ v2x = v2theta2 * v2x;
+ v2cx = v2ones + v2x;
+ v2sx = v2theta + v2theta * v2x;
+ v2cx[0] = v2sx[0];
+ __v4sf v4sx = _mm_cvtpd_ps (v2cx);
+ *sinx = v4sx[0];
+ *cosx = v4sx[1];
+ }
+ else
+ {
+ /* Handle some special cases. */
+ if (ix)
+ *sinx = theta - (theta * SMALL);
+ else
+ *sinx = theta;
+ *cosx = 1.0 - abstheta;
+ }
+ }
+ else /* |x| >= Pi/4. */
+ {
+ unsigned int signbit = xi >> 31;
+ if (ix < 0x40e231d6) /* |x| < 9*Pi/4. */
+ {
+ /* There are cases where FE_UPWARD rounding mode can
+ produce a result of abstheta * inv_PI_4 == 9,
+ where abstheta < 9pi/4, so the domain for
+ pio2_table must go to 5 (9 / 2 + 1). */
+ unsigned int n = (abstheta * inv_PI_4) + 1;
+ theta = abstheta - pio2_table[n / 2];
+ reduced_sincos (theta, n, signbit, sinx, cosx);
+ }
+ else if (ix < 0x7f800000)
+ {
+ if (ix < 0x4b000000) /* |x| < 2^23. */
+ {
+ unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
+ double x = n / 2;
+ theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
+ /* Argument reduction needed. */
+ reduced_sincos (theta, n, signbit, sinx, cosx);
+ }
+ else /* |x| >= 2^23. */
+ {
+ x = fabsf (x);
+ int exponent
+ = (ix >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS;
+ exponent += 3;
+ exponent /= 28;
+ double a = invpio4_table[exponent] * x;
+ double b = invpio4_table[exponent + 1] * x;
+ double c = invpio4_table[exponent + 2] * x;
+ double d = invpio4_table[exponent + 3] * x;
+ uint64_t l = a;
+ l &= ~0x7;
+ a -= l;
+ double e = a + b;
+ l = e;
+ e = a - l;
+ if (l & 1)
+ {
+ e -= 1.0;
+ e += b;
+ e += c;
+ e += d;
+ e *= M_PI_4;
+ reduced_sincos (e, l + 1, signbit, sinx, cosx);
+ }
+ else
+ {
+ e += b;
+ e += c;
+ e += d;
+ if (e <= 1.0)
+ {
+ e *= M_PI_4;
+ reduced_sincos (e, l + 1, signbit, sinx, cosx);
+ }
+ else
+ {
+ l++;
+ e -= 2.0;
+ e *= M_PI_4;
+ reduced_sincos (e, l + 1, signbit, sinx, cosx);
+ }
+ }
+ }
+ }
+ else
+ {
+ if (ix == 0x7f800000)
+ __set_errno (EDOM);
+ /* sin/cos(Inf or NaN) is NaN. */
+ *sinx = *cosx = x - x;
+ }
+ }
+}
+
+#ifndef SINCOSF
+libm_alias_float (__sincos, sincos)
+#endif
diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S b/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S
new file mode 100644
index 0000000000..51d012bb12
--- /dev/null
+++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S
@@ -0,0 +1,2 @@
+#define __sincosf __sincosf_sse2
+#include <sysdeps/x86_64/fpu/s_sincosf.S>
diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf.c b/sysdeps/x86_64/fpu/multiarch/s_sincosf.c
new file mode 100644
index 0000000000..6cb4295558
--- /dev/null
+++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf.c
@@ -0,0 +1,28 @@
+/* Multiple versions of sincosf.
+ Copyright (C) 2017 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <libm-alias-float.h>
+
+extern void __redirect_sincosf (float, float *, float *);
+
+#define SYMBOL_NAME sincosf
+#include "ifunc-fma.h"
+
+libc_ifunc_redirected (__redirect_sincosf, __sincosf, IFUNC_SELECTOR ());
+
+libm_alias_float (__sincos, sincos)
diff --git a/sysdeps/x86_64/fpu/s_sincosf.S b/sysdeps/x86_64/fpu/s_sincosf.S
index f00cd0eb3c..2086e8ca5c 100644
--- a/sysdeps/x86_64/fpu/s_sincosf.S
+++ b/sysdeps/x86_64/fpu/s_sincosf.S
@@ -561,4 +561,6 @@ L(SP_ONE):
.type L(SP_ONE), @object
ASM_SIZE_DIRECTIVE(L(SP_ONE))
+#ifndef __sincosf
libm_alias_float (__sincos, sincos)
+#endif
--
2.14.3