From patchwork Mon Jan 8 16:03:10 2018 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: "H.J. Lu" X-Patchwork-Id: 25271 Received: (qmail 130106 invoked by alias); 8 Jan 2018 16:03:16 -0000 Mailing-List: contact libc-alpha-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-alpha-owner@sourceware.org Delivered-To: mailing list libc-alpha@sourceware.org Received: (qmail 130096 invoked by uid 89); 8 Jan 2018 16:03:16 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-23.8 required=5.0 tests=AWL, BAYES_00, FREEMAIL_FROM, GIT_PATCH_0, GIT_PATCH_1, GIT_PATCH_2, GIT_PATCH_3, KAM_SHORT, KAM_STOCKGEN, RCVD_IN_DNSWL_NONE, SPF_PASS autolearn=ham version=3.3.2 spammy= X-HELO: mail-ot0-f196.google.com X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:mime-version:in-reply-to:references:from:date :message-id:subject:to; bh=l6YQ3eydYI9afj6aLhHO9G4iEUOAkOpgSnEFBOYKIFE=; b=cHbONlezdX7vItaAM2adNZJ6D8eQfmH5ywEvUe6qHqY+1k14D4GYTEGZjgs/I427xK 0fwrv9C5b0ylO5GqDk8C8C06pqiXu1VQiNiwMBreSEboP3nQpRiPmLCRmRoCBGmX0+XE +KZuWIB/QjqCBYq9vQ6v+25PX3UZSmbHykef+1O2k16DN8qUBXv9k1Q2xoG0uu7qK3Zd 8LvBHFx6aoo4wPxTvIOxVU3e6bRyd9zEEMuS3CJ4MecxNtG7vaJQWK7CjV4k265oZHK+ 1cMJJIrqGgG+CVJpmj5PHHYeYKpj7poe4FY3dWSQNOcCDXuGG/iyzELexgVagRW1YUUu OvMw== X-Gm-Message-State: AKwxyteY+0SkQWuZ7NpVkych5Khomo2l/3Wwu9XIdl8N8FtlCSqs9wJc Wa7TZGU4A2lUVCA8s+CqipKUAyQGV1/5Wol1Ay6QWQ== X-Google-Smtp-Source: ACJfBot9k0JVf8zOSs4s2Yr+96OJTsgqA1ZmvcdNf04UDbKbF4BFeJdGA6WpH+gfFTFWUbUsGH8cegrYrKBCVoRjMM0= X-Received: by 10.157.50.162 with SMTP id u31mr6726570otb.372.1515427391152; Mon, 08 Jan 2018 08:03:11 -0800 (PST) MIME-Version: 1.0 In-Reply-To: References: <20171220220054.GA16094@intel.com> <20171221003115.GA2338@intel.com> From: "H.J. Lu" Date: Mon, 8 Jan 2018 08:03:10 -0800 Message-ID: Subject: Re: [PATCH] x86-64: Add sincosf with vector FMA To: GNU C Library On Mon, Jan 8, 2018 at 4:43 AM, H.J. Lu wrote: > On Wed, Dec 20, 2017 at 4:31 PM, H.J. Lu 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 . 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% From 21b5bca0c20b55525adb1fe8373a61e87e441b5f Mon Sep 17 00:00:00 2001 From: "H.J. Lu" 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 + . */ + +#include +#include +#include +#include +#include +#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 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 + . */ + +#include + +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