This is the mail archive of the
glibc-cvs@sourceware.org
mailing list for the glibc project.
GNU C Library master sources branch master updated. glibc-2.28.9000-492-gba4b8fa
- From: hjl at sourceware dot org
- To: glibc-cvs at sourceware dot org
- Date: 26 Dec 2018 15:37:14 -0000
- Subject: GNU C Library master sources branch master updated. glibc-2.28.9000-492-gba4b8fa
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU C Library master sources".
The branch, master has been updated
via ba4b8fab20d52e4c3c52fcff805c77219a9e8b5e (commit)
via 9412979a4384c183eb50dbc6fe40aaf9904b1dc7 (commit)
via 8700a7851bccce63335f937f930382de58b8a249 (commit)
from 57b3ff8e1a4ca8118c5d9a47584e6e0ce6f7002f (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=ba4b8fab20d52e4c3c52fcff805c77219a9e8b5e
commit ba4b8fab20d52e4c3c52fcff805c77219a9e8b5e
Author: H.J. Lu <hjl.tools@gmail.com>
Date: Wed Dec 26 06:58:15 2018 -0800
x86-64: Remove s_sincosf-sse2.S
The current s_sincosf.c is faster than s_sincosf-sse2.S. On Broadwell
with FMA disabled, bench-sincosf shows:
Before After Improvement
max 154.032 114.517 34%
min 6.25 5.609 11%
mean 14.8728 12.8589 15%
* sysdeps/x86_64/fpu/s_sincosf.S: Removed.
* sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S: Likewise.
* sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.c: New file.
diff --git a/ChangeLog b/ChangeLog
index 42f6523..cb959e3 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,11 @@
2018-12-26 H.J. Lu <hongjiu.lu@intel.com>
+ * sysdeps/x86_64/fpu/s_sincosf.S: Removed.
+ * sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.c: New file.
+
+2018-12-26 H.J. Lu <hongjiu.lu@intel.com>
+
* sysdeps/x86_64/fpu/libm-test-ulps: Regenerated.
2018-12-26 H.J. Lu <hongjiu.lu@intel.com>
diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S b/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S
deleted file mode 100644
index 51d012b..0000000
--- a/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S
+++ /dev/null
@@ -1,2 +0,0 @@
-#define __sincosf __sincosf_sse2
-#include <sysdeps/x86_64/fpu/s_sincosf.S>
diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.c b/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.c
new file mode 100644
index 0000000..3c6d293
--- /dev/null
+++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.c
@@ -0,0 +1,2 @@
+#define SINCOSF __sincosf_sse2
+#include <sysdeps/ieee754/flt-32/s_sincosf.c>
diff --git a/sysdeps/x86_64/fpu/s_sincosf.S b/sysdeps/x86_64/fpu/s_sincosf.S
deleted file mode 100644
index 2086e8c..0000000
--- a/sysdeps/x86_64/fpu/s_sincosf.S
+++ /dev/null
@@ -1,566 +0,0 @@
-/* Optimized sincosf function.
- Copyright (C) 2012-2018 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 <sysdep.h>
-#include <errno.h>
-#include <libm-alias-float.h>
-
-/* Short algorithm description:
- *
- * 1) if |x|==0: sin(x)=x,
- * cos(x)=1.
- * 2) if |x|<2^-27: sin(x)=x-x*DP_SMALL, raising underflow only when needed,
- * cos(x)=1-|x|.
- * 3) if |x|<2^-5 : sin(x)=x+x*x^2*DP_SIN2_0+x^5*DP_SIN2_1,
- * cos(x)=1+1*x^2*DP_COS2_0+x^5*DP_COS2_1
- * 4) if |x|< Pi/4: sin(x)=x+x*x^2*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))),
- * cos(x)=1+1*x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))).
- * 5) if |x| < 9*Pi/4:
- * 5.1) Range reduction:
- * k=trunc(|x|/(Pi/4)), j=(k+1)&0x0e, n=k+1, t=|x|-j*Pi/4.
- * 5.2) Reconstruction:
- * sign_sin = sign(x) * (-1.0)^(( n >>2)&1)
- * sign_cos = (-1.0)^(((n+2)>>2)&1)
- * poly_sin = ((((S4*t^2 + S3)*t^2 + S2)*t^2 + S1)*t^2 + S0)*t^2*t+t
- * poly_cos = ((((C4*t^2 + C3)*t^2 + C2)*t^2 + C1)*t^2 + C0)*t^2*s+s
- * if(n&2 != 0) {
- * using cos(t) and sin(t) polynomials for |t|<Pi/4, results are
- * cos(x) = poly_sin * sign_cos
- * sin(x) = poly_cos * sign_sin
- * } else {
- * sin(x) = poly_sin * sign_sin
- * cos(x) = poly_cos * sign_cos
- * }
- * 6) if |x| < 2^23, large args:
- * 6.1) Range reduction:
- * k=trunc(|x|/(Pi/4)), j=(k+1)&0xfffffffe, n=k+1, t=|x|-j*Pi/4
- * 6.2) Reconstruction same as (5.2).
- * 7) if |x| >= 2^23, very large args:
- * 7.1) Range reduction:
- * k=trunc(|x|/(Pi/4)), j=(k+1)&0xfffffffe, n=k+1, t=|x|-j*Pi/4.
- * 7.2) Reconstruction same as (5.2).
- * 8) if x is Inf, return x-x, and set errno=EDOM.
- * 9) if x is NaN, return x-x.
- *
- * Special cases:
- * sin/cos(+-0) = +-0/1 not raising inexact/underflow,
- * sin/cos(subnormal) raises inexact/underflow,
- * sin/cos(min_normalized) raises inexact/underflow,
- * sin/cos(normalized) raises inexact,
- * sin/cos(Inf) = NaN, raises invalid, sets errno to EDOM,
- * sin/cos(NaN) = NaN.
- */
-
-# define ARG_SIN_PTR %rdi
-# define ARG_COS_PTR %rsi
-
- .text
-ENTRY(__sincosf)
- /* Input: %xmm0 contains single precision argument x */
- /* %rdi points to sin result */
- /* %rsi points to cos result */
-
- movd %xmm0, %eax /* Bits of x */
- movaps %xmm0, %xmm7 /* Copy of x */
- cvtss2sd %xmm0, %xmm0 /* DP x */
- movss L(SP_ABS_MASK)(%rip), %xmm3
- movl %eax, %r8d /* Copy of x bits */
- andl $0x7fffffff, %eax /* |x| */
-
- cmpl $0x3f490fdb, %eax /* |x|<Pi/4 ? */
- jb L(arg_less_pio4)
-
- /* Here if |x|>=Pi/4 */
- andps %xmm7, %xmm3 /* SP |x| */
- andpd L(DP_ABS_MASK)(%rip),%xmm0 /* DP |x| */
- movss L(SP_INVPIO4)(%rip), %xmm2 /* SP 1/(Pi/4) */
-
- cmpl $0x40e231d6, %eax /* |x|<9*Pi/4 ? */
- jae L(large_args)
-
- /* Here if Pi/4<=|x|<9*Pi/4 */
- mulss %xmm3, %xmm2 /* SP |x|/(Pi/4) */
- movl %r8d, %ecx /* Load x */
- cvttss2si %xmm2, %eax /* k, number of Pi/4 in x */
- lea L(PIO4J)(%rip), %r9
- shrl $29, %ecx /* (sign of x) << 2 */
- addl $1, %eax /* k+1 */
- movl $0x0e, %edx
- andl %eax, %edx /* j = (k+1)&0x0e */
- subsd (%r9,%rdx,8), %xmm0 /* t = |x| - j * Pi/4 */
-
-L(reconstruction):
- /* Input: %eax=n, %xmm0=t, %ecx=sign(x) */
-
- movaps %xmm0, %xmm4 /* t */
- movhpd L(DP_ONES)(%rip), %xmm4 /* 1|t */
- mulsd %xmm0, %xmm0 /* y=t^2 */
- movl $2, %edx
- unpcklpd %xmm0, %xmm0 /* y|y */
- addl %eax, %edx /* k+2 */
- movaps %xmm0, %xmm1 /* y|y */
- mulpd %xmm0, %xmm0 /* z=t^4|z=t^4 */
-
- movaps L(DP_SC4)(%rip), %xmm2 /* S4 */
- mulpd %xmm0, %xmm2 /* z*S4 */
- movaps L(DP_SC3)(%rip), %xmm3 /* S3 */
- mulpd %xmm0, %xmm3 /* z*S3 */
- xorl %eax, %ecx /* (sign_x ^ (k>>2))<<2 */
- addpd L(DP_SC2)(%rip), %xmm2 /* S2+z*S4 */
- mulpd %xmm0, %xmm2 /* z*(S2+z*S4) */
- shrl $2, %edx /* (k+2)>>2 */
- addpd L(DP_SC1)(%rip), %xmm3 /* S1+z*S3 */
- mulpd %xmm0, %xmm3 /* z*(S1+z*S3) */
- shrl $2, %ecx /* sign_x ^ k>>2 */
- addpd L(DP_SC0)(%rip), %xmm2 /* S0+z*(S2+z*S4) */
- andl $1, %edx /* sign_cos = ((k+2)>>2)&1 */
- mulpd %xmm1, %xmm2 /* y*(S0+z*(S2+z*S4)) */
- andl $1, %ecx /* sign_sin = sign_x ^ ((k>>2)&1) */
- addpd %xmm2, %xmm3 /* y*(S0+y*(S1+y*(S2+y*(S3+y*S4)))) */
- lea L(DP_ONES)(%rip), %r9
- mulpd %xmm4, %xmm3 /*t*y*(S0+y*(S1+y*(S2+y*(S3+y*S4))))*/
- testl $2, %eax /* n&2 != 0 ? */
- addpd %xmm4, %xmm3 /*t+t*y*(S0+y*(S1+y*(S2+y*(S3+y*S4))*/
- jnz L(sin_result_sin_poly)
-
-/*L(sin_result_cos_poly):*/
- /*
- * Here if
- * cos(x) = poly_sin * sign_cos
- * sin(x) = poly_cos * sign_sin
- */
- movsd (%r9,%rcx,8), %xmm4 /* 0|sign_sin */
- movhpd (%r9,%rdx,8), %xmm4 /* sign_cos|sign_sin */
- mulpd %xmm4, %xmm3 /* result_cos|result_sin */
- cvtpd2ps %xmm3, %xmm0 /* SP results */
- movss %xmm0, (ARG_SIN_PTR) /* store sin(x) from xmm0[0] */
- shufps $1, %xmm0, %xmm0 /* move cos(x) to xmm0[0] */
- movss %xmm0, (ARG_COS_PTR) /* store cos(x) */
- ret
-
- .p2align 4
-L(sin_result_sin_poly):
- /*
- * Here if
- * sin(x) = poly_sin * sign_sin
- * cos(x) = poly_cos * sign_cos
- */
- movsd (%r9,%rdx,8), %xmm4 /* 0|sign_cos */
- movhpd (%r9,%rcx,8), %xmm4 /* sign_sin|sign_cos */
- mulpd %xmm4, %xmm3 /* result_sin|result_cos */
- cvtpd2ps %xmm3, %xmm0 /* SP results */
- movss %xmm0, (ARG_COS_PTR) /* store cos(x) from xmm0[0] */
- shufps $1, %xmm0, %xmm0 /* move sin(x) to xmm0[0] */
- movss %xmm0, (ARG_SIN_PTR) /* store sin(x) */
- ret
-
- .p2align 4
-L(large_args):
- /* Here if |x|>=9*Pi/4 */
- cmpl $0x7f800000, %eax /* x is Inf or NaN ? */
- jae L(arg_inf_or_nan)
-
- /* Here if finite |x|>=9*Pi/4 */
- cmpl $0x4b000000, %eax /* |x|<2^23 ? */
- jae L(very_large_args)
-
- /* Here if 9*Pi/4<=|x|<2^23 */
- movsd L(DP_INVPIO4)(%rip), %xmm1 /* 1/(Pi/4) */
- mulsd %xmm0, %xmm1 /* |x|/(Pi/4) */
- cvttsd2si %xmm1, %eax /* k=trunc(|x|/(Pi/4)) */
- addl $1, %eax /* k+1 */
- movl %eax, %edx
- andl $0xfffffffe, %edx /* j=(k+1)&0xfffffffe */
- cvtsi2sdl %edx, %xmm4 /* DP j */
- movl %r8d, %ecx /* Load x */
- movsd L(DP_PIO4HI)(%rip), %xmm2 /* -PIO4HI = high part of -Pi/4 */
- shrl $29, %ecx /* (sign of x) << 2 */
- mulsd %xmm4, %xmm2 /* -j*PIO4HI */
- movsd L(DP_PIO4LO)(%rip), %xmm3 /* -PIO4LO = low part of -Pi/4 */
- addsd %xmm2, %xmm0 /* |x| - j*PIO4HI */
- mulsd %xmm3, %xmm4 /* j*PIO4LO */
- addsd %xmm4, %xmm0 /* t = |x| - j*PIO4HI - j*PIO4LO */
- jmp L(reconstruction)
-
- .p2align 4
-L(very_large_args):
- /* Here if finite |x|>=2^23 */
-
- /* bitpos = (ix>>23) - BIAS_32 + 59; */
- shrl $23, %eax /* eb = biased exponent of x */
- subl $68, %eax /* bitpos=eb-0x7f+59, where 0x7f */
- /*is exponent bias */
- movl $28, %ecx /* %cl=28 */
- movl %eax, %edx /* bitpos copy */
-
- /* j = bitpos/28; */
- div %cl /* j in register %al=%ax/%cl */
- movapd %xmm0, %xmm3 /* |x| */
- andl $0xff, %eax /* clear unneeded remainder from %ah*/
-
- imull $28, %eax, %ecx /* j*28 */
- lea L(_FPI)(%rip), %r9
- movsd L(DP_HI_MASK)(%rip), %xmm4 /* DP_HI_MASK */
- movapd %xmm0, %xmm5 /* |x| */
- mulsd -16(%r9,%rax,8), %xmm3 /* tmp3 = FPI[j-2]*|x| */
- movapd %xmm0, %xmm1 /* |x| */
- mulsd -8(%r9,%rax,8), %xmm5 /* tmp2 = FPI[j-1]*|x| */
- mulsd (%r9,%rax,8), %xmm0 /* tmp0 = FPI[j]*|x| */
- addl $19, %ecx /* j*28+19 */
- mulsd 8(%r9,%rax,8), %xmm1 /* tmp1 = FPI[j+1]*|x| */
- cmpl %ecx, %edx /* bitpos>=j*28+19 ? */
- jl L(very_large_skip1)
-
- /* Here if bitpos>=j*28+19 */
- andpd %xmm3, %xmm4 /* HI(tmp3) */
- subsd %xmm4, %xmm3 /* tmp3 = tmp3 - HI(tmp3) */
-L(very_large_skip1):
-
- movsd L(DP_2POW52)(%rip), %xmm6
- movapd %xmm5, %xmm2 /* tmp2 copy */
- addsd %xmm3, %xmm5 /* tmp5 = tmp3 + tmp2 */
- movl $1, %edx
- addsd %xmm5, %xmm6 /* tmp6 = tmp5 + 2^52 */
- movsd 8+L(DP_2POW52)(%rip), %xmm4
- movd %xmm6, %eax /* k = I64_LO(tmp6); */
- addsd %xmm6, %xmm4 /* tmp4 = tmp6 - 2^52 */
- movl %r8d, %ecx /* Load x */
- comisd %xmm5, %xmm4 /* tmp4 > tmp5 ? */
- jbe L(very_large_skip2)
-
- /* Here if tmp4 > tmp5 */
- subl $1, %eax /* k-- */
- addsd 8+L(DP_ONES)(%rip), %xmm4 /* tmp4 -= 1.0 */
-L(very_large_skip2):
-
- andl %eax, %edx /* k&1 */
- lea L(DP_ZERONE)(%rip), %r9
- subsd %xmm4, %xmm3 /* tmp3 -= tmp4 */
- addsd (%r9,%rdx,8), %xmm3 /* t = DP_ZERONE[k&1] + tmp3 */
- addsd %xmm2, %xmm3 /* t += tmp2 */
- shrl $29, %ecx /* (sign of x) << 2 */
- addsd %xmm3, %xmm0 /* t += tmp0 */
- addl $1, %eax /* n=k+1 */
- addsd %xmm1, %xmm0 /* t += tmp1 */
- mulsd L(DP_PIO4)(%rip), %xmm0 /* t *= PI04 */
-
- jmp L(reconstruction) /* end of very_large_args peth */
-
- .p2align 4
-L(arg_less_pio4):
- /* Here if |x|<Pi/4 */
- cmpl $0x3d000000, %eax /* |x|<2^-5 ? */
- jl L(arg_less_2pn5)
-
- /* Here if 2^-5<=|x|<Pi/4 */
- movaps %xmm0, %xmm3 /* DP x */
- movhpd L(DP_ONES)(%rip), %xmm3 /* DP 1|x */
- mulsd %xmm0, %xmm0 /* DP y=x^2 */
- unpcklpd %xmm0, %xmm0 /* DP y|y */
- movaps %xmm0, %xmm1 /* y|y */
- mulpd %xmm0, %xmm0 /* z=x^4|z=x^4 */
-
- movapd L(DP_SC4)(%rip), %xmm4 /* S4 */
- mulpd %xmm0, %xmm4 /* z*S4 */
- movapd L(DP_SC3)(%rip), %xmm5 /* S3 */
- mulpd %xmm0, %xmm5 /* z*S3 */
- addpd L(DP_SC2)(%rip), %xmm4 /* S2+z*S4 */
- mulpd %xmm0, %xmm4 /* z*(S2+z*S4) */
- addpd L(DP_SC1)(%rip), %xmm5 /* S1+z*S3 */
- mulpd %xmm0, %xmm5 /* z*(S1+z*S3) */
- addpd L(DP_SC0)(%rip), %xmm4 /* S0+z*(S2+z*S4) */
- mulpd %xmm1, %xmm4 /* y*(S0+z*(S2+z*S4)) */
- mulpd %xmm3, %xmm5 /* x*z*(S1+z*S3) */
- mulpd %xmm3, %xmm4 /* x*y*(S0+z*(S2+z*S4)) */
- addpd %xmm5, %xmm4 /*x*y*(S0+y*(S1+y*(S2+y*(S3+y*S4))))*/
- addpd %xmm4, %xmm3 /*x+x*y*(S0+y*(S1+y*(S2+y*(S3+y*S4))*/
- cvtpd2ps %xmm3, %xmm0 /* SP results */
- movss %xmm0, (ARG_SIN_PTR) /* store sin(x) from xmm0[0] */
- shufps $1, %xmm0, %xmm0 /* move cos(x) to xmm0[0] */
- movss %xmm0, (ARG_COS_PTR) /* store cos(x) */
- ret
-
- .p2align 4
-L(arg_less_2pn5):
- /* Here if |x|<2^-5 */
- cmpl $0x32000000, %eax /* |x|<2^-27 ? */
- jl L(arg_less_2pn27)
-
- /* Here if 2^-27<=|x|<2^-5 */
- movaps %xmm0, %xmm1 /* DP x */
- movhpd L(DP_ONES)(%rip), %xmm1 /* DP 1|x */
- mulsd %xmm0, %xmm0 /* DP x^2 */
- unpcklpd %xmm0, %xmm0 /* DP x^2|x^2 */
-
- movaps L(DP_SINCOS2_1)(%rip), %xmm3 /* DP DP_SIN2_1 */
- mulpd %xmm0, %xmm3 /* DP x^2*DP_SIN2_1 */
- addpd L(DP_SINCOS2_0)(%rip), %xmm3 /* DP DP_SIN2_0+x^2*DP_SIN2_1 */
- mulpd %xmm0, %xmm3 /* DP x^2*DP_SIN2_0+x^4*DP_SIN2_1 */
- mulpd %xmm1, %xmm3 /* DP x^3*DP_SIN2_0+x^5*DP_SIN2_1 */
- addpd %xmm1, %xmm3 /* DP x+x^3*DP_SIN2_0+x^5*DP_SIN2_1 */
- cvtpd2ps %xmm3, %xmm0 /* SP results */
- movss %xmm0, (ARG_SIN_PTR) /* store sin(x) from xmm0[0] */
- shufps $1, %xmm0, %xmm0 /* move cos(x) to xmm0[0] */
- movss %xmm0, (ARG_COS_PTR) /* store cos(x) */
- ret
-
- .p2align 4
-L(arg_less_2pn27):
- cmpl $0, %eax /* x=0 ? */
- je L(arg_zero) /* in case x=0 return sin(+-0)==+-0 */
- /* Here if |x|<2^-27 */
- /*
- * Special cases here:
- * sin(subnormal) raises inexact/underflow
- * sin(min_normalized) raises inexact/underflow
- * sin(normalized) raises inexact
- * cos(here)=1-|x| (raising inexact)
- */
- movaps %xmm0, %xmm3 /* DP x */
- mulsd L(DP_SMALL)(%rip), %xmm0/* DP x*DP_SMALL */
- subsd %xmm0, %xmm3 /* DP sin result is x-x*DP_SMALL */
- andps L(SP_ABS_MASK)(%rip), %xmm7/* SP |x| */
- cvtsd2ss %xmm3, %xmm0 /* sin(x) */
- movss L(SP_ONE)(%rip), %xmm1 /* SP 1.0 */
- movss %xmm0, (ARG_SIN_PTR) /* sin(x) store */
- subss %xmm7, %xmm1 /* cos(x) */
- movss %xmm1, (ARG_COS_PTR) /* cos(x) store */
- ret
-
- .p2align 4
-L(arg_zero):
- movss L(SP_ONE)(%rip), %xmm0 /* 1.0 */
- movss %xmm7, (ARG_SIN_PTR) /* sin(+-0)==x */
- movss %xmm0, (ARG_COS_PTR) /* cos(+-0)==1 */
- ret
-
- .p2align 4
-L(arg_inf_or_nan):
- /* Here if |x| is Inf or NAN */
- jne L(skip_errno_setting) /* in case of x is NaN */
-
- /* Align stack to 16 bytes. */
- subq $8, %rsp
- cfi_adjust_cfa_offset (8)
- /* Here if x is Inf. Set errno to EDOM. */
- call JUMPTARGET(__errno_location)
- addq $8, %rsp
- cfi_adjust_cfa_offset (-8)
-
- movl $EDOM, (%rax)
-
- .p2align 4
-L(skip_errno_setting):
- /* Here if |x| is Inf or NAN. Continued. */
- subss %xmm7, %xmm7 /* x-x, result is NaN */
- movss %xmm7, (ARG_SIN_PTR)
- movss %xmm7, (ARG_COS_PTR)
- ret
-END(__sincosf)
-
- .section .rodata, "a"
- .p2align 3
-L(PIO4J): /* Table of j*Pi/4, for j=0,1,..,10 */
- .long 0x00000000,0x00000000
- .long 0x54442d18,0x3fe921fb
- .long 0x54442d18,0x3ff921fb
- .long 0x7f3321d2,0x4002d97c
- .long 0x54442d18,0x400921fb
- .long 0x2955385e,0x400f6a7a
- .long 0x7f3321d2,0x4012d97c
- .long 0xe9bba775,0x4015fdbb
- .long 0x54442d18,0x401921fb
- .long 0xbeccb2bb,0x401c463a
- .long 0x2955385e,0x401f6a7a
- .type L(PIO4J), @object
- ASM_SIZE_DIRECTIVE(L(PIO4J))
-
- .p2align 3
-L(_FPI): /* 4/Pi broken into sum of positive DP values */
- .long 0x00000000,0x00000000
- .long 0x6c000000,0x3ff45f30
- .long 0x2a000000,0x3e3c9c88
- .long 0xa8000000,0x3c54fe13
- .long 0xd0000000,0x3aaf47d4
- .long 0x6c000000,0x38fbb81b
- .long 0xe0000000,0x3714acc9
- .long 0x7c000000,0x3560e410
- .long 0x56000000,0x33bca2c7
- .long 0xac000000,0x31fbd778
- .long 0xe0000000,0x300b7246
- .long 0xe8000000,0x2e5d2126
- .long 0x48000000,0x2c970032
- .long 0xe8000000,0x2ad77504
- .long 0xe0000000,0x290921cf
- .long 0xb0000000,0x274deb1c
- .long 0xe0000000,0x25829a73
- .long 0xbe000000,0x23fd1046
- .long 0x10000000,0x2224baed
- .long 0x8e000000,0x20709d33
- .long 0x80000000,0x1e535a2f
- .long 0x64000000,0x1cef904e
- .long 0x30000000,0x1b0d6398
- .long 0x24000000,0x1964ce7d
- .long 0x16000000,0x17b908bf
- .type L(_FPI), @object
- ASM_SIZE_DIRECTIVE(L(_FPI))
-
-/* Coefficients of polynomials for */
-/* sin(x)~=x+x*x^2*(DP_SIN2_0+x^2*DP_SIN2_1) in low DP part, */
-/* cos(x)~=1+1*x^2*(DP_COS2_0+x^2*DP_COS2_1) in high DP part, */
-/* for |x|<2^-5. */
- .p2align 4
-L(DP_SINCOS2_0):
- .long 0x5543d49d,0xbfc55555
- .long 0xff5cc6fd,0xbfdfffff
- .type L(DP_SINCOS2_0), @object
- ASM_SIZE_DIRECTIVE(L(DP_SINCOS2_0))
-
- .p2align 4
-L(DP_SINCOS2_1):
- .long 0x75cec8c5,0x3f8110f4
- .long 0xb178dac5,0x3fa55514
- .type L(DP_SINCOS2_1), @object
- ASM_SIZE_DIRECTIVE(L(DP_SINCOS2_1))
-
-
- .p2align 3
-L(DP_ZERONE):
- .long 0x00000000,0x00000000 /* 0.0 */
- .long 0x00000000,0xbff00000 /* 1.0 */
- .type L(DP_ZERONE), @object
- ASM_SIZE_DIRECTIVE(L(DP_ZERONE))
-
- .p2align 3
-L(DP_ONES):
- .long 0x00000000,0x3ff00000 /* +1.0 */
- .long 0x00000000,0xbff00000 /* -1.0 */
- .type L(DP_ONES), @object
- ASM_SIZE_DIRECTIVE(L(DP_ONES))
-
-/* Coefficients of polynomials for */
-/* sin(t)~=t+t*t^2*(S0+t^2*(S1+t^2*(S2+t^2*(S3+t^2*S4)))) in low DP part, */
-/* cos(t)~=1+1*t^2*(C0+t^2*(C1+t^2*(C2+t^2*(C3+t^2*C4)))) in high DP part, */
-/* for |t|<Pi/4. */
- .p2align 4
-L(DP_SC4):
- .long 0x1674b58a,0xbe5a947e
- .long 0xdd8844d7,0xbe923c97
- .type L(DP_SC4), @object
- ASM_SIZE_DIRECTIVE(L(DP_SC4))
-
- .p2align 4
-L(DP_SC3):
- .long 0x64e6b5b4,0x3ec71d72
- .long 0x9ac43cc0,0x3efa00eb
- .type L(DP_SC3), @object
- ASM_SIZE_DIRECTIVE(L(DP_SC3))
-
- .p2align 4
-L(DP_SC2):
- .long 0x8b4bd1f9,0xbf2a019f
- .long 0x348b6874,0xbf56c16b
- .type L(DP_SC2), @object
- ASM_SIZE_DIRECTIVE(L(DP_SC2))
-
- .p2align 4
-L(DP_SC1):
- .long 0x10c2688b,0x3f811111
- .long 0x545c50c7,0x3fa55555
- .type L(DP_SC1), @object
- ASM_SIZE_DIRECTIVE(L(DP_SC1))
-
- .p2align 4
-L(DP_SC0):
- .long 0x55551cd9,0xbfc55555
- .long 0xfffe98ae,0xbfdfffff
- .type L(DP_SC0), @object
- ASM_SIZE_DIRECTIVE(L(DP_SC0))
-
- .p2align 3
-L(DP_SMALL):
- .long 0x00000000,0x3cd00000 /* 2^(-50) */
- .type L(DP_SMALL), @object
- ASM_SIZE_DIRECTIVE(L(DP_SMALL))
-
- .p2align 3
-L(DP_PIO4):
- .long 0x54442d18,0x3fe921fb /* Pi/4 */
- .type L(DP_PIO4), @object
- ASM_SIZE_DIRECTIVE(L(DP_PIO4))
-
- .p2align 3
-L(DP_2POW52):
- .long 0x00000000,0x43300000 /* +2^52 */
- .long 0x00000000,0xc3300000 /* -2^52 */
- .type L(DP_2POW52), @object
- ASM_SIZE_DIRECTIVE(L(DP_2POW52))
-
- .p2align 3
-L(DP_INVPIO4):
- .long 0x6dc9c883,0x3ff45f30 /* 4/Pi */
- .type L(DP_INVPIO4), @object
- ASM_SIZE_DIRECTIVE(L(DP_INVPIO4))
-
- .p2align 3
-L(DP_PIO4HI):
- .long 0x54000000,0xbfe921fb /* High part of Pi/4 */
- .type L(DP_PIO4HI), @object
- ASM_SIZE_DIRECTIVE(L(DP_PIO4HI))
-
- .p2align 3
-L(DP_PIO4LO):
- .long 0x11A62633,0xbe010b46 /* Low part of Pi/4 */
- .type L(DP_PIO4LO), @object
- ASM_SIZE_DIRECTIVE(L(DP_PIO4LO))
-
- .p2align 2
-L(SP_INVPIO4):
- .long 0x3fa2f983 /* 4/Pi */
- .type L(SP_INVPIO4), @object
- ASM_SIZE_DIRECTIVE(L(SP_INVPIO4))
-
- .p2align 4
-L(DP_ABS_MASK): /* Mask for getting DP absolute value */
- .long 0xffffffff,0x7fffffff
- .long 0xffffffff,0x7fffffff
- .type L(DP_ABS_MASK), @object
- ASM_SIZE_DIRECTIVE(L(DP_ABS_MASK))
-
- .p2align 3
-L(DP_HI_MASK): /* Mask for getting high 21 bits of DP value */
- .long 0x00000000,0xffffffff
- .type L(DP_HI_MASK), @object
- ASM_SIZE_DIRECTIVE(L(DP_HI_MASK))
-
- .p2align 4
-L(SP_ABS_MASK): /* Mask for getting SP absolute value */
- .long 0x7fffffff,0x7fffffff
- .long 0x7fffffff,0x7fffffff
- .type L(SP_ABS_MASK), @object
- ASM_SIZE_DIRECTIVE(L(SP_ABS_MASK))
-
- .p2align 2
-L(SP_ONE):
- .long 0x3f800000 /* 1.0 */
- .type L(SP_ONE), @object
- ASM_SIZE_DIRECTIVE(L(SP_ONE))
-
-#ifndef __sincosf
-libm_alias_float (__sincos, sincos)
-#endif
http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=9412979a4384c183eb50dbc6fe40aaf9904b1dc7
commit 9412979a4384c183eb50dbc6fe40aaf9904b1dc7
Author: H.J. Lu <hjl.tools@gmail.com>
Date: Wed Dec 26 06:57:17 2018 -0800
Regenerate sysdeps/x86_64/fpu/libm-test-ulps
* sysdeps/x86_64/fpu/libm-test-ulps: Regenerated.
diff --git a/ChangeLog b/ChangeLog
index 6767bb7..42f6523 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,9 @@
2018-12-26 H.J. Lu <hongjiu.lu@intel.com>
+ * sysdeps/x86_64/fpu/libm-test-ulps: Regenerated.
+
+2018-12-26 H.J. Lu <hongjiu.lu@intel.com>
+
* sysdeps/ieee754/flt-32/s_sincosf.h: Include <sincosf_poly.h>.
(sincos_t, sincosf_poly, sinf_poly): Moved to ...
* sysdeps/ieee754/flt-32/sincosf_poly.h: Here. New file.
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 493ad87..ee5bcbd 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -2603,24 +2603,30 @@ ldouble: 1
Function: "sincos_downward":
double: 1
+float: 1
float128: 3
idouble: 1
+ifloat: 1
ifloat128: 3
ildouble: 3
ldouble: 3
Function: "sincos_towardzero":
double: 1
+float: 1
float128: 2
idouble: 1
+ifloat: 1
ifloat128: 2
ildouble: 2
ldouble: 2
Function: "sincos_upward":
double: 1
+float: 1
float128: 3
idouble: 1
+ifloat: 1
ifloat128: 3
ildouble: 3
ldouble: 3
http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=8700a7851bccce63335f937f930382de58b8a249
commit 8700a7851bccce63335f937f930382de58b8a249
Author: H.J. Lu <hjl.tools@gmail.com>
Date: Wed Dec 26 06:56:04 2018 -0800
x86-64: Vectorize sincosf_poly and update s_sincosf-fma.c
Add <sincosf_poly.h> and include it in s_sincosf.h to allow vectorized
sincosf_poly. Add x86 sincosf_poly.h to vectorize sincosf_poly. On
Broadwell, bench-sincosf shows:
Before After Improvement
max 160.273 114.198 40%
min 6.25 5.625 11%
mean 13.0325 10.6462 22%
Vectorized sincosf_poly shows
Before After Improvement
max 138.653 114.198 21%
min 5.004 5.625 -11%
mean 11.5934 10.6462 9%
Tested on x86-64 and i686 as well as with build-many-glibcs.py.
* sysdeps/ieee754/flt-32/s_sincosf.h: Include <sincosf_poly.h>.
(sincos_t, sincosf_poly, sinf_poly): Moved to ...
* sysdeps/ieee754/flt-32/sincosf_poly.h: Here. New file.
* sysdeps/x86/fpu/s_sincosf_data.c: New file.
* sysdeps/x86/fpu/sincosf_poly.h: Likewise.
* sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: Just include
<sysdeps/ieee754/flt-32/s_sincosf.c>.
diff --git a/ChangeLog b/ChangeLog
index 881808d..6767bb7 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+2018-12-26 H.J. Lu <hongjiu.lu@intel.com>
+
+ * sysdeps/ieee754/flt-32/s_sincosf.h: Include <sincosf_poly.h>.
+ (sincos_t, sincosf_poly, sinf_poly): Moved to ...
+ * sysdeps/ieee754/flt-32/sincosf_poly.h: Here. New file.
+ * sysdeps/x86/fpu/s_sincosf_data.c: New file.
+ * sysdeps/x86/fpu/sincosf_poly.h: Likewise.
+ * sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: Just include
+ <sysdeps/ieee754/flt-32/s_sincosf.c>.
+
2018-12-21 Joseph Myers <joseph@codesourcery.com>
[BZ #24023]
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h
index 1dcb04f..54a9da5 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.h
+++ b/sysdeps/ieee754/flt-32/s_sincosf.h
@@ -19,22 +19,13 @@
#include <stdint.h>
#include <math.h>
#include "math_config.h"
+#include <sincosf_poly.h>
/* 2PI * 2^-64. */
static const double pi63 = 0x1.921FB54442D18p-62;
/* PI / 4. */
static const double pio4 = 0x1.921FB54442D18p-1;
-/* The constants and polynomials for sine and cosine. */
-typedef struct
-{
- double sign[4]; /* Sign of sine in quadrants 0..3. */
- double hpi_inv; /* 2 / PI ( * 2^24 if !TOINT_INTRINSICS). */
- double hpi; /* PI / 2. */
- double c0, c1, c2, c3, c4; /* Cosine polynomial. */
- double s1, s2, s3; /* Sine polynomial. */
-} sincos_t;
-
/* Polynomial data (the cosine polynomial is negated in the 2nd entry). */
extern const sincos_t __sincosf_table[2] attribute_hidden;
@@ -48,66 +39,6 @@ abstop12 (float x)
return (asuint (x) >> 20) & 0x7ff;
}
-/* Compute the sine and cosine of inputs X and X2 (X squared), using the
- polynomial P and store the results in SINP and COSP. N is the quadrant,
- if odd the cosine and sine polynomials are swapped. */
-static inline void
-sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp,
- float *cosp)
-{
- double x3, x4, x5, x6, s, c, c1, c2, s1;
-
- x4 = x2 * x2;
- x3 = x2 * x;
- c2 = p->c3 + x2 * p->c4;
- s1 = p->s2 + x2 * p->s3;
-
- /* Swap sin/cos result based on quadrant. */
- float *tmp = (n & 1 ? cosp : sinp);
- cosp = (n & 1 ? sinp : cosp);
- sinp = tmp;
-
- c1 = p->c0 + x2 * p->c1;
- x5 = x3 * x2;
- x6 = x4 * x2;
-
- s = x + x3 * p->s1;
- c = c1 + x4 * p->c2;
-
- *sinp = s + x5 * s1;
- *cosp = c + x6 * c2;
-}
-
-/* Return the sine of inputs X and X2 (X squared) using the polynomial P.
- N is the quadrant, and if odd the cosine polynomial is used. */
-static inline float
-sinf_poly (double x, double x2, const sincos_t *p, int n)
-{
- double x3, x4, x6, x7, s, c, c1, c2, s1;
-
- if ((n & 1) == 0)
- {
- x3 = x * x2;
- s1 = p->s2 + x2 * p->s3;
-
- x7 = x3 * x2;
- s = x + x3 * p->s1;
-
- return s + x7 * s1;
- }
- else
- {
- x4 = x2 * x2;
- c2 = p->c3 + x2 * p->c4;
- c1 = p->c0 + x2 * p->c1;
-
- x6 = x4 * x2;
- c = c1 + x4 * p->c2;
-
- return c + x6 * c2;
- }
-}
-
/* Fast range reduction using single multiply-subtract. Return the modulo of
X as a value between -PI/4 and PI/4 and store the quadrant in NP.
The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double
diff --git a/sysdeps/ieee754/flt-32/sincosf_poly.h b/sysdeps/ieee754/flt-32/sincosf_poly.h
new file mode 100644
index 0000000..f49c75a
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/sincosf_poly.h
@@ -0,0 +1,87 @@
+/* Used by sinf, cosf and sincosf functions.
+ Copyright (C) 2018 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/>. */
+
+/* The constants and polynomials for sine and cosine. */
+typedef struct
+{
+ double sign[4]; /* Sign of sine in quadrants 0..3. */
+ double hpi_inv; /* 2 / PI ( * 2^24 if !TOINT_INTRINSICS). */
+ double hpi; /* PI / 2. */
+ double c0, c1, c2, c3, c4; /* Cosine polynomial. */
+ double s1, s2, s3; /* Sine polynomial. */
+} sincos_t;
+
+/* Compute the sine and cosine of inputs X and X2 (X squared), using the
+ polynomial P and store the results in SINP and COSP. N is the quadrant,
+ if odd the cosine and sine polynomials are swapped. */
+static inline void
+sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp,
+ float *cosp)
+{
+ double x3, x4, x5, x6, s, c, c1, c2, s1;
+
+ x4 = x2 * x2;
+ x3 = x2 * x;
+ c2 = p->c3 + x2 * p->c4;
+ s1 = p->s2 + x2 * p->s3;
+
+ /* Swap sin/cos result based on quadrant. */
+ float *tmp = (n & 1 ? cosp : sinp);
+ cosp = (n & 1 ? sinp : cosp);
+ sinp = tmp;
+
+ c1 = p->c0 + x2 * p->c1;
+ x5 = x3 * x2;
+ x6 = x4 * x2;
+
+ s = x + x3 * p->s1;
+ c = c1 + x4 * p->c2;
+
+ *sinp = s + x5 * s1;
+ *cosp = c + x6 * c2;
+}
+
+/* Return the sine of inputs X and X2 (X squared) using the polynomial P.
+ N is the quadrant, and if odd the cosine polynomial is used. */
+static inline float
+sinf_poly (double x, double x2, const sincos_t *p, int n)
+{
+ double x3, x4, x6, x7, s, c, c1, c2, s1;
+
+ if ((n & 1) == 0)
+ {
+ x3 = x * x2;
+ s1 = p->s2 + x2 * p->s3;
+
+ x7 = x3 * x2;
+ s = x + x3 * p->s1;
+
+ return s + x7 * s1;
+ }
+ else
+ {
+ x4 = x2 * x2;
+ c2 = p->c3 + x2 * p->c4;
+ c1 = p->c0 + x2 * p->c1;
+
+ x6 = x4 * x2;
+ c = c1 + x4 * p->c2;
+
+ return c + x6 * c2;
+ }
+}
diff --git a/sysdeps/x86/fpu/s_sincosf_data.c b/sysdeps/x86/fpu/s_sincosf_data.c
new file mode 100644
index 0000000..b219644
--- /dev/null
+++ b/sysdeps/x86/fpu/s_sincosf_data.c
@@ -0,0 +1,68 @@
+/* Compute sine and cosine of argument.
+ Copyright (C) 2018 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 <stdint.h>
+#include <math.h>
+#include <sysdeps/ieee754/flt-32/math_config.h>
+#include <s_sincosf.h>
+
+/* The constants and polynomials for sine and cosine. The 2nd entry
+ computes -cos (x) rather than cos (x) to get negation for free. */
+const sincos_t __sincosf_table[2] =
+{
+ {
+ { 1.0, -1.0, -1.0, 1.0 },
+#if TOINT_INTRINSICS
+ 0x1.45F306DC9C883p-1,
+#else
+ 0x1.45F306DC9C883p+23,
+#endif
+ 0x1.921FB54442D18p0,
+ 0x1p0,
+ -0x1.ffffffd0c621cp-2,
+ { -0x1.555545995a603p-3, 0x1.55553e1068f19p-5 },
+ { 0x1.1107605230bc4p-7, -0x1.6c087e89a359dp-10 },
+ { -0x1.994eb3774cf24p-13, 0x1.99343027bf8c3p-16 }
+ },
+ {
+ { 1.0, -1.0, -1.0, 1.0 },
+#if TOINT_INTRINSICS
+ 0x1.45F306DC9C883p-1,
+#else
+ 0x1.45F306DC9C883p+23,
+#endif
+ 0x1.921FB54442D18p0,
+ -0x1p0,
+ 0x1.ffffffd0c621cp-2,
+ { -0x1.555545995a603p-3, -0x1.55553e1068f19p-5 },
+ { 0x1.1107605230bc4p-7, 0x1.6c087e89a359dp-10 },
+ { -0x1.994eb3774cf24p-13, -0x1.99343027bf8c3p-16 }
+ }
+};
+
+/* Table with 4/PI to 192 bit precision. To avoid unaligned accesses
+ only 8 new bits are added per entry, making the table 4 times larger. */
+const uint32_t __inv_pio4[24] =
+{
+ 0xa2, 0xa2f9, 0xa2f983, 0xa2f9836e,
+ 0xf9836e4e, 0x836e4e44, 0x6e4e4415, 0x4e441529,
+ 0x441529fc, 0x1529fc27, 0x29fc2757, 0xfc2757d1,
+ 0x2757d1f5, 0x57d1f534, 0xd1f534dd, 0xf534ddc0,
+ 0x34ddc0db, 0xddc0db62, 0xc0db6295, 0xdb629599,
+ 0x6295993c, 0x95993c43, 0x993c4390, 0x3c439041
+};
diff --git a/sysdeps/x86/fpu/sincosf_poly.h b/sysdeps/x86/fpu/sincosf_poly.h
new file mode 100644
index 0000000..396ba89
--- /dev/null
+++ b/sysdeps/x86/fpu/sincosf_poly.h
@@ -0,0 +1,111 @@
+/* Used by sinf, cosf and sincosf functions. X86-64 version.
+ Copyright (C) 2018 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/>. */
+
+typedef double v2df_t __attribute__ ((vector_size (2 * sizeof (double))));
+
+#ifdef __SSE2_MATH__
+typedef float v4sf_t __attribute__ ((vector_size (4 * sizeof (float))));
+
+static inline void
+v2df_to_sf (v2df_t v2df, float *f0p, float *f1p)
+{
+ v4sf_t v4sf = __builtin_ia32_cvtpd2ps (v2df);
+ *f0p = v4sf[0];
+ *f1p = v4sf[1];
+}
+#else
+static inline void
+v2df_to_sf (v2df_t v2df, float *f0p, float *f1p)
+{
+ *f0p = (float) v2df[0];
+ *f1p = (float) v2df[1];
+}
+#endif
+
+/* The constants and polynomials for sine and cosine. */
+typedef struct
+{
+ double sign[4]; /* Sign of sine in quadrants 0..3. */
+ double hpi_inv; /* 2 / PI ( * 2^24 if !TOINT_INTRINSICS). */
+ double hpi; /* PI / 2. */
+ /* Cosine polynomial: c0, c1, c2, c3, c4.
+ Sine polynomial: s1, s2, s3. */
+ double c0, c1;
+ v2df_t s1c2, s2c3, s3c4;
+} sincos_t;
+
+/* Compute the sine and cosine of inputs X and X2 (X squared), using the
+ polynomial P and store the results in SINP and COSP. N is the quadrant,
+ if odd the cosine and sine polynomials are swapped. */
+static inline void
+sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp,
+ float *cosp)
+{
+ v2df_t vx2x2 = { x2, x2 };
+ v2df_t vxx2 = { x, x2 };
+ v2df_t vx3x4, vs1c2;
+
+ vx3x4 = vx2x2 * vxx2;
+ vs1c2 = p->s2c3 + vx2x2 * p->s3c4;
+
+ /* Swap sin/cos result based on quadrant. */
+ if (n & 1)
+ {
+ float *tmp = cosp;
+ cosp = sinp;
+ sinp = tmp;
+ }
+
+ double c1 = p->c0 + x2 * p->c1;
+ v2df_t vxc1 = { x, c1 };
+ v2df_t vx5x6 = vx3x4 * vx2x2;
+
+ v2df_t vsincos = vxc1 + vx3x4 * p->s1c2;
+ vsincos = vsincos + vx5x6 * vs1c2;
+ v2df_to_sf (vsincos, sinp, cosp);
+}
+
+/* Return the sine of inputs X and X2 (X squared) using the polynomial P.
+ N is the quadrant, and if odd the cosine polynomial is used. */
+static inline float
+sinf_poly (double x, double x2, const sincos_t *p, int n)
+{
+ double x3, x4, x6, x7, s, c, c1, c2, s1;
+
+ if ((n & 1) == 0)
+ {
+ x3 = x * x2;
+ s1 = p->s2c3[0] + x2 * p->s3c4[0];
+
+ x7 = x3 * x2;
+ s = x + x3 * p->s1c2[0];
+
+ return s + x7 * s1;
+ }
+ else
+ {
+ x4 = x2 * x2;
+ c2 = p->s2c3[1] + x2 * p->s3c4[1];
+ c1 = p->c0 + x2 * p->c1;
+
+ x6 = x4 * x2;
+ c = c1 + x4 * p->s1c2[1];
+
+ return c + x6 * c2;
+ }
+}
diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c
index 0b80c4f..253dab1 100644
--- a/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c
+++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c
@@ -1,271 +1,2 @@
-/* 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>
-
#define SINCOSF __sincosf_fma
-
-#ifndef SINCOSF
-# define SINCOSF_FUNC __sincosf
-#else
-# define SINCOSF_FUNC SINCOSF
-#endif
-
-/* PI/2 with 98 bits of accuracy. */
-static const double PI_2_hi = 0x1.921fb544p+0;
-static const double PI_2_lo = 0x1.0b4611a626332p-34;
-
-static const double SMALL = 0x1p-50; /* 2^-50. */
-static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */
-
-#define FLOAT_EXPONENT_SHIFT 23
-#define FLOAT_EXPONENT_BIAS 127
-
-static const double pio2_table[] = {
- 0 * M_PI_2,
- 1 * M_PI_2,
- 2 * M_PI_2,
- 3 * M_PI_2,
- 4 * M_PI_2,
- 5 * M_PI_2
-};
-
-static const double invpio4_table[] = {
- 0x0p+0,
- 0x1.45f306cp+0,
- 0x1.c9c882ap-28,
- 0x1.4fe13a8p-58,
- 0x1.f47d4dp-85,
- 0x1.bb81b6cp-112,
- 0x1.4acc9ep-142,
- 0x1.0e4107cp-169
-};
-
-static const double ones[] = { 1.0, -1.0 };
-
-/* 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
+#include <sysdeps/ieee754/flt-32/s_sincosf.c>
-----------------------------------------------------------------------
Summary of changes:
ChangeLog | 20 +
sysdeps/ieee754/flt-32/s_sincosf.h | 71 +---
sysdeps/ieee754/flt-32/sincosf_poly.h | 87 ++++
sysdeps/x86/fpu/s_sincosf_data.c | 68 +++
sysdeps/x86/fpu/sincosf_poly.h | 111 +++++
sysdeps/x86_64/fpu/libm-test-ulps | 6 +
sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c | 271 +------------
sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S | 2 -
sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.c | 2 +
sysdeps/x86_64/fpu/s_sincosf.S | 566 -------------------------
10 files changed, 296 insertions(+), 908 deletions(-)
create mode 100644 sysdeps/ieee754/flt-32/sincosf_poly.h
create mode 100644 sysdeps/x86/fpu/s_sincosf_data.c
create mode 100644 sysdeps/x86/fpu/sincosf_poly.h
delete mode 100644 sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S
create mode 100644 sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.c
delete mode 100644 sysdeps/x86_64/fpu/s_sincosf.S
hooks/post-receive
--
GNU C Library master sources