diff mbox series

[v3,12/12] math: Evaluate the sincosf sine and cosine as two-lane vectors

Message ID 20260923215118.2608685-13-adhemerval.zanella@linaro.org (mailing list archive)
State Under Review
Delegated to: Wilco Dijkstra
Headers
Series Add sinf/cosf/sincosf CORE-MATH implementations |

Checks

Context Check Description
redhat-pt-bot/TryBot-apply_patch success Patch applied to master at the time it was sent
linaro-tcwg-bot/tcwg_glibc_build--master-arm success Build passed
redhat-pt-bot/TryBot-32bit success Build for i686
linaro-tcwg-bot/tcwg_glibc_build--master-aarch64 success Build passed
linaro-tcwg-bot/tcwg_glibc_check--master-arm success Test passed
linaro-tcwg-bot/tcwg_glibc_check--master-aarch64 success Test passed

Commit Message

Adhemerval Zanella Sept. 23, 2026, 9:50 p.m. UTC 
  Compute the sine and the cosine of the sincosf fast path together with
the GCC generic vector extension.  It keeps the same computation from
the original CORE-MATH, so the results are unchanged.

The polynomial coefficients are interleaved, and the sine and cosine
tables become a single table of { cos, -sin } pairs, from which the fast
path loads { c0, -s0 } and { -s0, -c0 } as aligned vectors.  The sinf
and cosf read the same tables, which shrink from 1664 to 704 bytes.

Targets with double-precision vector units evaluate both polynomials
and the final combination in one pass.  On aarch64 with gcc 16 the
fast path goes from 58 to 48 instructions.

Improvement over the previous commit in percent, reciprocal-throughput /
latency (same builds and machines as the CORE-MATH import commits):

                   random       moderate           huge        blender
x86_64        -1.7 / -1.7    +0.0 / +4.2    -1.0 / -0.8    -1.1 / -0.9
x86_64v2      -1.3 / -1.1    +3.2 / +9.5    +0.1 / +0.0    -0.6 / -0.4
x86_64v3      -0.8 / -1.6    +8.9 / +8.5    -0.3 / -0.2    -0.8 / -1.5
aarch64     +11.4 / +11.5    +1.0 / +1.1    -0.4 / -0.5  +11.4 / +11.5
power10       +3.1 / +3.6    +0.3 / -0.4    +2.9 / +3.0    +1.5 / +1.8
loongarch     +1.0 / +1.0    -0.3 / -0.1    +1.0 / +1.1    +1.0 / +1.0

aarch64 gains about 11% on the random and blender workloads.

sinf and cosf, which now read the same tables (geomean over their
workloads):

                     sinf           cosf
x86_64        +1.8 / -0.7    -2.9 / -0.3
x86_64v2      +7.0 / +0.4    -3.8 / -1.0
x86_64v3      +6.4 / +0.3    -1.7 / +0.8
aarch64       +0.5 / +0.4    +0.8 / +1.3
power10       +4.9 / +2.1    +4.9 / +4.8
loongarch     +1.9 / +0.1    +2.1 / -1.4

Checked exhaustively (all inputs, all rounding modes) against the
correctly rounded results on x86_64, i686 (SSE2), aarch64, powerpc64le,
powerpc64, loongarch64, riscv64, and s390x.
---
 sysdeps/ieee754/flt-32/s_cosf.c         | 12 ++--
 sysdeps/ieee754/flt-32/s_sincosf.c      | 34 +++++++----
 sysdeps/ieee754/flt-32/s_sincosf_data.c | 81 ++++++++++++++-----------
 sysdeps/ieee754/flt-32/s_sincosf_data.h | 19 +++---
 sysdeps/ieee754/flt-32/s_sinf.c         | 16 +++--
 5 files changed, 96 insertions(+), 66 deletions(-)
diff mbox series

Patch

diff --git a/sysdeps/ieee754/flt-32/s_cosf.c b/sysdeps/ieee754/flt-32/s_cosf.c
index 2ca4d91694e..d625815f47a 100644
--- a/sysdeps/ieee754/flt-32/s_cosf.c
+++ b/sysdeps/ieee754/flt-32/s_cosf.c
@@ -59,9 +59,9 @@  as_cosf_big (float x)
   int ia;
   double z = RBIG_SINCOSF (t, &ia);
   double z2 = z * z, z4 = z2 * z2;
-  double aa = (A[0] + z2 * A[1]) + z4 * (A[2] + z2 * A[3]);
-  double bb = (B[0] + z2 * B[1]) + z4 * (B[2] + z2 * B[3]);
-  double s0 = TB_COSF[(ia + 8u) & 31], c0 = TB_COSF[ia & 31];
+  double aa = (A (0) + z2 * A (1)) + z4 * (A (2) + z2 * A (3));
+  double bb = (B (0) + z2 * B (1)) + z4 * (B (2) + z2 * B (3));
+  double c0 = CS[ia & 31][0], s0 = CS[ia & 31][1];
   double r = c0 + z * (aa * s0 - bb * (z * c0));
   uint32_t tr = asuint64 (r);
   uint64_t tail = (tr + 6) & (~UINT64_C(0) >> 36);
@@ -109,9 +109,9 @@  COSF_FUNC (float x)
       z = rltl (z0, &ia);
     }
   double z2 = z * z, z4 = z2 * z2;
-  double aa = (A[0] + z2 * A[1]) + z4 * (A[2] + z2 * A[3]);
-  double bb = (B[0] + z2 * B[1]) + z4 * (B[2] + z2 * B[3]);
-  double c0 = TB_COSF[ia & 31], s0 = TB_COSF[(ia + 8) & 31];
+  double aa = (A (0) + z2 * A (1)) + z4 * (A (2) + z2 * A (3));
+  double bb = (B (0) + z2 * B (1)) + z4 * (B (2) + z2 * B (3));
+  double c0 = CS[ia & 31][0], s0 = CS[ia & 31][1];
   double r = c0 + aa * (z * s0) - bb * (z2 * c0);
   return r;
 }
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c
index b55d0c7675e..15c0bdeb4d5 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.c
+++ b/sysdeps/ieee754/flt-32/s_sincosf.c
@@ -76,10 +76,12 @@  as_sincosf_big (float x, float *sout, float *cout)
   int ia;
   double z = RBIG_SINCOSF (t, &ia);
   double z2 = z * z, z4 = z2 * z2;
-  double aa = (A[0] + z2 * A[1]) + z4 * (A[2] + z2 * A[3]);
-  double bb = (B[0] + z2 * B[1]) + z4 * (B[2] + z2 * B[3]);
+  double aa = (A (0) + z2 * A (1)) + z4 * (A (2) + z2 * A (3));
+  double bb = (B (0) + z2 * B (1)) + z4 * (B (2) + z2 * B (3));
   bb *= z;
-  double s0 = TB[ia & 31], c0 = TB[(ia + 8u) & 31];
+  /* CS[(ia + 24) & 31] = { sin(ia*pi/16), cos(ia*pi/16) }.  */
+  unsigned int k = (ia + 24u) & 31;
+  double s0 = CS[k][0], c0 = CS[k][1];
   double s = s0 + z * (aa * c0 - bb * s0);
   double c = c0 - z * (aa * s0 + bb * c0);
   *sout = s;
@@ -133,15 +135,23 @@  SINCOSF_FUNC (float x, float *sout, float *cout)
       z = rltl (z0, &ia);
     }
   double z2 = z * z, z4 = z2 * z2;
-  double aa = (A[0] + z2 * A[1]) + z4 * (A[2] + z2 * A[3]);
-  double bb = (B[0] + z2 * B[1]) + z4 * (B[2] + z2 * B[3]);
-  aa *= z;
-  bb *= z2;
-  double s0 = TB[ia & 31], c0 = TB[(ia + 8) & 31];
-  double rs = s0 + (aa * c0 - bb * s0);
-  double rc = c0 - (aa * s0 + bb * c0);
-  *sout = rs;
-  *cout = rc;
+  /* The sine and the cosine are evaluated together, lane 0 and lane 1
+     of the vectors, as
+       sin = s0 + (aa * c0 - bb * s0)
+       cos = c0 - (aa * s0 + bb * c0)
+     with { aa, bb } = q, m = CS[ia] = { c0, -s0 } and
+     w = CS[ia + 8] = { -s0, -c0 }: t = aa * m + bb * w and r = t - w.
+     Each lane performs the operations of the scalar formula, the cosine
+     up to exact negations, so the results do not depend on whether the
+     compiler uses vector instructions.  */
+  sincosf_v2df_t p = (AB[0] + z2 * AB[1]) + z4 * (AB[2] + z2 * AB[3]);
+  sincosf_v2df_t q = p * (sincosf_v2df_t) { z, z2 };
+  const sincosf_v2df_t *cs = CS + (ia & 31);
+  sincosf_v2df_t m = cs[0], w = cs[8];
+  sincosf_v2df_t t = q[0] * m + q[1] * w;
+  sincosf_v2df_t r = t - w;
+  *sout = r[0];
+  *cout = r[1];
 }
 
 #ifndef SINCOSF
diff --git a/sysdeps/ieee754/flt-32/s_sincosf_data.c b/sysdeps/ieee754/flt-32/s_sincosf_data.c
index d248db8c08f..ef2f8192917 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf_data.c
+++ b/sysdeps/ieee754/flt-32/s_sincosf_data.c
@@ -35,44 +35,57 @@  const uint64_t __sinf_ipi[] =
 
 const sincosf_tables_t __sincosf_tables =
 {
-  .b =
+  .cs =
   {
-     0x1.3bd3cc9be45dcp-6, -0x1.03c1f081b0833p-14,  0x1.55d3c6fc9ac1fp-24,
-    -0x1.e1d3ff281b40dp-35
+    { 0x1p+0, 0x0p+0 },
+    { 0x1.f6297cff75cbp-1, -0x1.8f8b83c69a60bp-3 },
+    { 0x1.d906bcf328d46p-1, -0x1.87de2a6aea963p-2 },
+    { 0x1.a9b66290ea1a3p-1, -0x1.1c73b39ae68c8p-1 },
+    { 0x1.6a09e667f3bcdp-1, -0x1.6a09e667f3bcdp-1 },
+    { 0x1.1c73b39ae68c8p-1, -0x1.a9b66290ea1a3p-1 },
+    { 0x1.87de2a6aea963p-2, -0x1.d906bcf328d46p-1 },
+    { 0x1.8f8b83c69a60bp-3, -0x1.f6297cff75cbp-1 },
+    { 0x0p+0, -0x1p+0 },
+    { -0x1.8f8b83c69a60bp-3, -0x1.f6297cff75cbp-1 },
+    { -0x1.87de2a6aea963p-2, -0x1.d906bcf328d46p-1 },
+    { -0x1.1c73b39ae68c8p-1, -0x1.a9b66290ea1a3p-1 },
+    { -0x1.6a09e667f3bcdp-1, -0x1.6a09e667f3bcdp-1 },
+    { -0x1.a9b66290ea1a3p-1, -0x1.1c73b39ae68c8p-1 },
+    { -0x1.d906bcf328d46p-1, -0x1.87de2a6aea963p-2 },
+    { -0x1.f6297cff75cbp-1, -0x1.8f8b83c69a60bp-3 },
+    { -0x1p+0, 0x0p+0 },
+    { -0x1.f6297cff75cbp-1, 0x1.8f8b83c69a60bp-3 },
+    { -0x1.d906bcf328d46p-1, 0x1.87de2a6aea963p-2 },
+    { -0x1.a9b66290ea1a3p-1, 0x1.1c73b39ae68c8p-1 },
+    { -0x1.6a09e667f3bcdp-1, 0x1.6a09e667f3bcdp-1 },
+    { -0x1.1c73b39ae68c8p-1, 0x1.a9b66290ea1a3p-1 },
+    { -0x1.87de2a6aea963p-2, 0x1.d906bcf328d46p-1 },
+    { -0x1.8f8b83c69a60bp-3, 0x1.f6297cff75cbp-1 },
+    { 0x0p+0, 0x1p+0 },
+    { 0x1.8f8b83c69a60bp-3, 0x1.f6297cff75cbp-1 },
+    { 0x1.87de2a6aea963p-2, 0x1.d906bcf328d46p-1 },
+    { 0x1.1c73b39ae68c8p-1, 0x1.a9b66290ea1a3p-1 },
+    { 0x1.6a09e667f3bcdp-1, 0x1.6a09e667f3bcdp-1 },
+    { 0x1.a9b66290ea1a3p-1, 0x1.1c73b39ae68c8p-1 },
+    { 0x1.d906bcf328d46p-1, 0x1.87de2a6aea963p-2 },
+    { 0x1.f6297cff75cbp-1, 0x1.8f8b83c69a60bp-3 },
+    { 0x1p+0, 0x0p+0 },
+    { 0x1.f6297cff75cbp-1, -0x1.8f8b83c69a60bp-3 },
+    { 0x1.d906bcf328d46p-1, -0x1.87de2a6aea963p-2 },
+    { 0x1.a9b66290ea1a3p-1, -0x1.1c73b39ae68c8p-1 },
+    { 0x1.6a09e667f3bcdp-1, -0x1.6a09e667f3bcdp-1 },
+    { 0x1.1c73b39ae68c8p-1, -0x1.a9b66290ea1a3p-1 },
+    { 0x1.87de2a6aea963p-2, -0x1.d906bcf328d46p-1 },
+    { 0x1.8f8b83c69a60bp-3, -0x1.f6297cff75cbp-1 },
   },
-  .a =
+  .ab =
   {
-     0x1.921fb54442d17p-3, -0x1.4abbce6256a39p-10,   0x1.466bc5a518c16p-19,
-    -0x1.32bdc61074ff6p-29
-  },
-  .tb =
-  {
-     0x0p+0,                0x1.8f8b83c69a60bp-3,  0x1.87de2a6aea963p-2,
-     0x1.1c73b39ae68c8p-1,  0x1.6a09e667f3bcdp-1,  0x1.a9b66290ea1a3p-1,
-     0x1.d906bcf328d46p-1,  0x1.f6297cff75cbp-1,   0x1p+0,
-     0x1.f6297cff75cbp-1,   0x1.d906bcf328d46p-1,  0x1.a9b66290ea1a3p-1,
-     0x1.6a09e667f3bcdp-1,  0x1.1c73b39ae68c8p-1,  0x1.87de2a6aea963p-2,
-     0x1.8f8b83c69a60bp-3,  0x0p+0,               -0x1.8f8b83c69a60bp-3,
-    -0x1.87de2a6aea963p-2, -0x1.1c73b39ae68c8p-1, -0x1.6a09e667f3bcdp-1,
-    -0x1.a9b66290ea1a3p-1, -0x1.d906bcf328d46p-1, -0x1.f6297cff75cbp-1,
-    -0x1p+0,               -0x1.f6297cff75cbp-1,  -0x1.d906bcf328d46p-1,
-    -0x1.a9b66290ea1a3p-1, -0x1.6a09e667f3bcdp-1, -0x1.1c73b39ae68c8p-1,
-    -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60bp-3
-  },
-  .tb_cosf =
-  {
-     0x1p+0,                0x1.f6297cff75cbp-1,   0x1.d906bcf328d46p-1,
-     0x1.a9b66290ea1a3p-1,  0x1.6a09e667f3bcdp-1,  0x1.1c73b39ae68c8p-1,
-     0x1.87de2a6aea963p-2,  0x1.8f8b83c69a60bp-3,  0x0p+0,
-    -0x1.8f8b83c69a60bp-3, -0x1.87de2a6aea963p-2, -0x1.1c73b39ae68c8p-1,
-    -0x1.6a09e667f3bcdp-1, -0x1.a9b66290ea1a3p-1, -0x1.d906bcf328d46p-1,
-    -0x1.f6297cff75cbp-1,  -0x1p+0,               -0x1.f6297cff75cbp-1,
-    -0x1.d906bcf328d46p-1, -0x1.a9b66290ea1a3p-1, -0x1.6a09e667f3bcdp-1,
-    -0x1.1c73b39ae68c8p-1, -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60bp-3,
-     0x0p+0,                0x1.8f8b83c69a60bp-3,  0x1.87de2a6aea963p-2,
-     0x1.1c73b39ae68c8p-1,  0x1.6a09e667f3bcdp-1,  0x1.a9b66290ea1a3p-1,
-     0x1.d906bcf328d46p-1,  0x1.f6297cff75cbp-1
+    { 0x1.921fb54442d17p-3, 0x1.3bd3cc9be45dcp-6 },
+    { -0x1.4abbce6256a39p-10, -0x1.03c1f081b0833p-14 },
+    { 0x1.466bc5a518c16p-19, 0x1.55d3c6fc9ac1fp-24 },
+    { -0x1.32bdc61074ff6p-29, -0x1.e1d3ff281b40dp-35 },
   },
+
 };
 
 const sincosf_database_t __sinf_st[] =
diff --git a/sysdeps/ieee754/flt-32/s_sincosf_data.h b/sysdeps/ieee754/flt-32/s_sincosf_data.h
index 7d9a7028f39..1bde95cbd66 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf_data.h
+++ b/sysdeps/ieee754/flt-32/s_sincosf_data.h
@@ -24,18 +24,21 @@ 
 extern const uint64_t __sinf_ipi[] attribute_hidden;
 #define IPI __sinf_ipi
 
+typedef double sincosf_v2df_t __attribute__ ((__vector_size__ (16)));
+
 typedef struct
 {
-  double b[4];
-  double a[4];
-  double tb[32];
-  double tb_cosf[32];
+  /* { cos(k*pi/16), -sin(k*pi/16) }, with the first 8 entries repeated
+     at the end so CS[(i & 31) + 8] needs no wrap around.  */
+  sincosf_v2df_t cs[40];
+  /* The coefficients { a[k], b[k] } of the two polynomials in z^2.  */
+  sincosf_v2df_t ab[4];
 } sincosf_tables_t;
 extern const sincosf_tables_t __sincosf_tables attribute_hidden;
-#define B __sincosf_tables.b
-#define A __sincosf_tables.a
-#define TB __sincosf_tables.tb
-#define TB_COSF __sincosf_tables.tb_cosf
+#define AB __sincosf_tables.ab
+#define A(k) AB[k][0]
+#define B(k) AB[k][1]
+#define CS __sincosf_tables.cs
 
 typedef struct
 {
diff --git a/sysdeps/ieee754/flt-32/s_sinf.c b/sysdeps/ieee754/flt-32/s_sinf.c
index b4f2b04903a..472546fb17f 100644
--- a/sysdeps/ieee754/flt-32/s_sinf.c
+++ b/sysdeps/ieee754/flt-32/s_sinf.c
@@ -68,9 +68,11 @@  as_sinf_big (float x)
   int ia;
   double z = RBIG_SINCOSF (t, &ia);
   double z2 = z * z, z4 = z2 * z2;
-  double aa = (A[0] + z2 * A[1]) + z4 * (A[2] + z2 * A[3]);
-  double bb = (B[0] + z2 * B[1]) + z4 * (B[2] + z2 * B[3]);
-  double s0 = TB[ia & 31], c0 = TB[(ia + 8u) & 31];
+  double aa = (A (0) + z2 * A (1)) + z4 * (A (2) + z2 * A (3));
+  double bb = (B (0) + z2 * B (1)) + z4 * (B (2) + z2 * B (3));
+  /* CS[(ia + 24) & 31] = { sin(ia*pi/16), cos(ia*pi/16) }.  */
+  unsigned int k = (ia + 24u) & 31;
+  double s0 = CS[k][0], c0 = CS[k][1];
   double r = s0 + z * (aa * c0 - bb * (z * s0));
   return r;
 }
@@ -119,9 +121,11 @@  SINF_FUNC (float x)
       z = rltl (z0, &ia);
     }
   double z2 = z * z, z4 = z2 * z2;
-  double aa = (A[0] + z2 * A[1]) + z4 * (A[2] + z2 * A[3]);
-  double bb = (B[0] + z2 * B[1]) + z4 * (B[2] + z2 * B[3]);
-  double s0 = TB[ia & 31], c0 = TB[(ia + 8) & 31];
+  double aa = (A (0) + z2 * A (1)) + z4 * (A (2) + z2 * A (3));
+  double bb = (B (0) + z2 * B (1)) + z4 * (B (2) + z2 * B (3));
+  /* CS[(ia + 24) & 31] = { sin(ia*pi/16), cos(ia*pi/16) }.  */
+  unsigned int k = (ia + 24u) & 31;
+  double s0 = CS[k][0], c0 = CS[k][1];
   double r = s0 + aa * (z * c0) - bb * (z2 * s0);
   return r;
 }