AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-21 20:40:05 +00:00
Code Issues Packages Projects Releases Wiki Activity

x86-64: Add sincosf with vector FMA

Browse Source 
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.
This commit is contained in:
H.J. Lu 2018-01-08 08:04:26 -08:00
parent 579396ee08
commit c70e4e9c9e
7 changed files with 288 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
ChangeLog
View File

@ -1,3 +1,14 @@
2018-01-08 H.J. Lu <hongjiu.lu@intel.com>
* 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.
2018-01-08 Florian Weimer <fweimer@redhat.com>
* nptl/tst-thread-exit-clobber.cc: New file.

4
NEWS
View File

@ -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.

5
sysdeps/x86_64/fpu/multiarch/Makefile
View File

@ -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 \

240
sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c Normal file
View File

@ -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

2
sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S Normal file
View File

@ -0,0 +1,2 @@
#define __sincosf __sincosf_sse2
#include <sysdeps/x86_64/fpu/s_sincosf.S>

28
sysdeps/x86_64/fpu/multiarch/s_sincosf.c Normal file
View File

@ -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)

2
sysdeps/x86_64/fpu/s_sincosf.S
View File

@ -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