[PATCH 7/8] math: Consolidate common definitions for cosf/sinf/tanf/sincosf

[Adhemerval Zanella <adhemerval.zanella@linaro.org>]

The rbig does not use floating-point instructions, so it can be on a
different TU and not be built multiple times for multiarch (since the
resulting code will be similar and not a hotspot).

Checked on x86_64-linux-gnu, i686-linux-gnu, aarch64-linux-gnu,
and powerpc64le-linux-gnu.
---
 math/Makefile                             |  1 +
 sysdeps/ieee754/flt-32/s_cosf.c           | 64 +---------------
 sysdeps/ieee754/flt-32/s_sincosf.c        | 71 +-----------------
 sysdeps/ieee754/flt-32/s_sincosf_common.c | 90 +++++++++++++++++++++++
 sysdeps/ieee754/flt-32/s_sincosf_common.h | 63 ++++++++++++++++
 sysdeps/ieee754/flt-32/s_sinf.c           | 71 +-----------------
 sysdeps/ieee754/flt-32/s_tanf.c           | 65 ++--------------
 7 files changed, 166 insertions(+), 259 deletions(-)
 create mode 100644 sysdeps/ieee754/flt-32/s_sincosf_common.c
 create mode 100644 sysdeps/ieee754/flt-32/s_sincosf_common.h

diff --git a/math/Makefile b/math/Makefile
index 7a9352c2cd..eaea4c77b9 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -389,6 +389,7 @@ type-float-routines := \
   s_asincosf_data \
   s_asincoshf_data \
   s_asincospif_data \
+  s_sincosf_common \
   s_sincosf_data \
   s_sincospif_data \
   # type-float-routines
diff --git a/sysdeps/ieee754/flt-32/s_cosf.c b/sysdeps/ieee754/flt-32/s_cosf.c
index 2bef3b3036..47ebc530c6 100644
--- a/sysdeps/ieee754/flt-32/s_cosf.c
+++ b/sysdeps/ieee754/flt-32/s_cosf.c
@@ -31,69 +31,9 @@ SOFTWARE.
 #include <libm-alias-float.h>
 #include "math_config.h"
 #include <math_uint128.h>
+#include <s_sincosf_common.h>
 #include <s_sincosf_data.h>
 
-static double __attribute__ ((noinline))
-rbig (uint32_t u, int *q)
-{
-  int e = (u >> 23) & 0xff, i;
-  uint64_t m = (u & (~0u >> 9)) | 1 << 23;
-  u128 p0 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[0]));
-  u128 p1 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[1]));
-  p1 = u128_add (p1, u128_rshift (p0, 64));
-  u128 p2 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[2]));
-  p2 = u128_add (p2, u128_rshift (p1, 64));
-  u128 p3 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[3]));
-  p3 = u128_add (p3, u128_rshift (p2, 64));
-  uint64_t p3h = u128_high (p3), p3l = u128_low (p3), p2l = u128_low (p2),
-	   p1l = u128_low (p1);
-  int64_t A;
-  int k = e - 124, s = k - 23;
-  /* in cr_sinf(), rbig() is called in the case 127+28 <= e < 0xff
-     thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
-  if (s < 64)
-    {
-      i = p3h << s | p3l >> (64 - s);
-      A = p3l << s | p2l >> (64 - s);
-    }
-  else if (s == 64)
-    {
-      i = p3l;
-      A = p2l;
-    }
-  else
-    { /* s > 64 */
-      i = p3l << (s - 64) | p2l >> (128 - s);
-      A = p2l << (s - 64) | p1l >> (128 - s);
-    }
-  int sgn = u;
-  sgn >>= 31;
-  int64_t sm = A >> 63;
-  i -= sm;
-  double z = (A ^ sgn) * 0x1p-64;
-  i = (i ^ sgn) - sgn;
-  *q = i;
-  return z;
-}
-
-static inline double
-rltl (float z, int *q)
-{
-  double x = z;
-  double idl = -0x1.b1bbead603d8bp-29 * x, idh = 0x1.45f306ep+2 * x,
-	 id = roundeven_finite (idh);
-  *q = asuint64 (0x1.8p52 + id);
-  return (idh - id) + idl;
-}
-
-static inline double
-rltl0 (double x, int *q)
-{
-  double idh = 0x1.45f306dc9c883p+2 * x, id = roundeven_finite (idh);
-  *q = asuint64 (0x1.8p52 + id);
-  return idh - id;
-}
-
 static float __attribute__ ((noinline))
 as_cosf_database (float x, double r)
 {
@@ -117,7 +57,7 @@ as_cosf_big (float x)
       return __math_invalidf (x);
     }
   int ia;
-  double z = rbig (t, &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]);
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c
index eb5d849174..9bb5b84a02 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.c
+++ b/sysdeps/ieee754/flt-32/s_sincosf.c
@@ -29,6 +29,7 @@ SOFTWARE.
 #include <libm-alias-float.h>
 #include "math_config.h"
 #include <math_uint128.h>
+#include <s_sincosf_common.h>
 #include <s_sincosf_data.h>
 
 #ifndef SECTION
@@ -41,74 +42,6 @@ SOFTWARE.
 # define SINCOSF_FUNC SINCOSF
 #endif
 
-static double __attribute__ ((noinline))
-rbig (uint32_t u, int *q)
-{
-  int e = (u >> 23) & 0xff, i;
-  uint64_t m = (u & (~0u >> 9)) | 1 << 23;
-  u128 p0 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[0]));
-  u128 p1 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[1]));
-  p1 = u128_add (p1, u128_rshift (p0, 64));
-  u128 p2 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[2]));
-  p2 = u128_add (p2, u128_rshift (p1, 64));
-  u128 p3 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[3]));
-  p3 = u128_add (p3, u128_rshift (p2, 64));
-  uint64_t p3h = u128_high (p3), p3l = u128_low (p3), p2l = u128_low (p2),
-	   p1l = u128_low (p1);
-  int64_t a;
-  int k = e - 124, s = k - 23;
-  /* in cr_sinf(), rbig() is called in the case 127+28 <= e < 0xff
-     thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
-  if (s < 64)
-    {
-      i = p3h << s | p3l >> (64 - s);
-      a = p3l << s | p2l >> (64 - s);
-    }
-  else if (s == 64)
-    {
-      i = p3l;
-      a = p2l;
-    }
-  else
-    { /* s > 64 */
-      i = p3l << (s - 64) | p2l >> (128 - s);
-      a = p2l << (s - 64) | p1l >> (128 - s);
-    }
-  int sgn = u;
-  sgn >>= 31;
-  int64_t sm = a >> 63;
-  i -= sm;
-  double z = (a ^ sgn) * 0x1p-64;
-  i = (i ^ sgn) - sgn;
-  *q = i;
-  return z;
-}
-
-static inline double
-rltl (float z, int *q)
-{
-  double x = z;
-  double idl = -0x1.b1bbead603d8bp-29 * x, idh = 0x1.45f306ep+2 * x,
-	 id = roundeven_finite (idh);
-  *q = asuint64 (0x1.8p52 + id);
-  return (idh - id) + idl;
-}
-
-static inline double
-rltl0 (double x, int *q)
-{
-  double idh = 0x1.45f306dc9c883p+2 * x, id = roundeven_finite (idh);
-  *q = asuint64 (0x1.8p52 + id);
-  return idh - id;
-}
-
-static inline float
-add_sign (float x, float rh, float rl)
-{
-  float sgn = copysignf (1.0f, x);
-  return sgn * rh + sgn * rl;
-}
-
 static void __attribute__ ((noinline))
 as_sincosf_database (float x, float *sout, float *cout)
 {
@@ -141,7 +74,7 @@ as_sincosf_big (float x, float *sout, float *cout)
       return;
     }
   int ia;
-  double z = rbig (t, &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]);
diff --git a/sysdeps/ieee754/flt-32/s_sincosf_common.c b/sysdeps/ieee754/flt-32/s_sincosf_common.c
new file mode 100644
index 0000000000..6d0f47d402
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/s_sincosf_common.c
@@ -0,0 +1,90 @@
+/* Common routines for sinf/cosf/tanf/sincosf.
+
+Copyright (c) 2022-2026 Alexei Sibidanov.
+
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/tan/tanf.c, revision 59d21d7).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+#include <math_uint128.h>
+#include <s_sincosf_common.h>
+#include <s_sincosf_data.h>
+
+static double
+rbig (uint32_t u, int *q, int e_bias)
+{
+  int e = (u >> 23) & 0xff, i;
+  uint64_t m = (u & (~0u >> 9)) | 1 << 23;
+  u128 p0 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[0]));
+  u128 p1 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[1]));
+  p1 = u128_add (p1, u128_rshift (p0, 64));
+  u128 p2 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[2]));
+  p2 = u128_add (p2, u128_rshift (p1, 64));
+  u128 p3 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[3]));
+  p3 = u128_add (p3, u128_rshift (p2, 64));
+  uint64_t p3h = u128_high (p3), p3l = u128_low (p3), p2l = u128_low (p2),
+	   p1l = u128_low (p1);
+  int64_t a;
+  int k = e - e_bias, s = k - 23;
+  /* in cr_sinf(), rbig() is called in the case 127+28 <= e < 0xff
+     thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
+  if (s < 64)
+    {
+      i = p3h << s | p3l >> (64 - s);
+      a = p3l << s | p2l >> (64 - s);
+    }
+  else if (s == 64)
+    {
+      i = p3l;
+      a = p2l;
+    }
+  else
+    { /* s > 64 */
+      i = p3l << (s - 64) | p2l >> (128 - s);
+      a = p2l << (s - 64) | p1l >> (128 - s);
+    }
+  int sgn = u;
+  sgn >>= 31;
+  int64_t sm = a >> 63;
+  i -= sm;
+  double z = (a ^ sgn) * 0x1p-64;
+  i = (i ^ sgn) - sgn;
+  *q = i;
+  return z;
+}
+
+double
+__sincosf_rbig (uint32_t u, int *q)
+{
+  /* in cr_sinf(), rbig() is called in the case 127+28 <= e < 0xff
+     thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
+  return rbig (u, q, 124);
+}
+
+/* argument reduction
+   same as rltl, but for |x| >= 2^28  */
+double
+__tanf_rbig (uint32_t u, int *q)
+{
+  /* in ctanf(), rbig() is called in the case 127+28 <= e < 0xff
+     thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
+  return rbig (u, q, 127);
+}
diff --git a/sysdeps/ieee754/flt-32/s_sincosf_common.h b/sysdeps/ieee754/flt-32/s_sincosf_common.h
new file mode 100644
index 0000000000..a15a0afe80
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/s_sincosf_common.h
@@ -0,0 +1,63 @@
+/* Common routines for sinf/cosf/tanf/sincosf.
+
+Copyright (c) 2022-2026 Alexei Sibidanov.
+
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/tan/tanf.c, revision 59d21d7).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+#ifndef _S_SINCOSF_COMMON_H
+#define _S_SINCOSF_COMMON_H
+
+#include <stdint.h>
+#include "math_config.h"
+
+static inline double
+rltl0 (double x, int *q)
+{
+  double idh = 0x1.45f306dc9c883p+2 * x, id = roundeven_finite (idh);
+  *q = asuint64 (0x1.8p52 + id);
+  return idh - id;
+}
+
+static inline double
+rltl (float z, int *q)
+{
+  double x = z;
+  double idl = -0x1.b1bbead603d8bp-29 * x, idh = 0x1.45f306ep+2 * x,
+	 id = roundeven_finite (idh);
+  *q = asuint64 (0x1.8p52 + id);
+  return (idh - id) + idl;
+}
+
+static inline float
+add_sign (float x, float rh, float rl)
+{
+  float sgn = copysignf (1.0f, x);
+  return sgn * rh + sgn * rl;
+}
+
+extern double __tanf_rbig (uint32_t u, int *q) attribute_hidden;
+#define RBIG_TANF __tanf_rbig
+extern double __sincosf_rbig (uint32_t u, int *q) attribute_hidden;
+#define RBIG_SINCOSF __sincosf_rbig
+
+#endif
diff --git a/sysdeps/ieee754/flt-32/s_sinf.c b/sysdeps/ieee754/flt-32/s_sinf.c
index 45d8a62a44..795a9131ba 100644
--- a/sysdeps/ieee754/flt-32/s_sinf.c
+++ b/sysdeps/ieee754/flt-32/s_sinf.c
@@ -31,6 +31,7 @@ SOFTWARE.
 #include <libm-alias-float.h>
 #include "math_config.h"
 #include <math_uint128.h>
+#include <s_sincosf_common.h>
 #include <s_sincosf_data.h>
 
 #ifndef SECTION
@@ -43,74 +44,6 @@ SOFTWARE.
 # define SINF_FUNC SINF
 #endif
 
-static double __attribute__ ((noinline))
-rbig (uint32_t u, int *q)
-{
-  int e = (u >> 23) & 0xff, i;
-  uint64_t m = (u & (~0u >> 9)) | 1 << 23;
-  u128 p0 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[0]));
-  u128 p1 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[1]));
-  p1 = u128_add (p1, u128_rshift (p0, 64));
-  u128 p2 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[2]));
-  p2 = u128_add (p2, u128_rshift (p1, 64));
-  u128 p3 = u128_mul (u128_from_u64 (m), u128_from_u64 (IPI[3]));
-  p3 = u128_add (p3, u128_rshift (p2, 64));
-  uint64_t p3h = u128_high (p3), p3l = u128_low (p3), p2l = u128_low (p2),
-	   p1l = u128_low (p1);
-  int64_t a;
-  int k = e - 124, s = k - 23;
-  /* in cr_sinf(), rbig() is called in the case 127+28 <= e < 0xff
-     thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
-  if (s < 64)
-    {
-      i = p3h << s | p3l >> (64 - s);
-      a = p3l << s | p2l >> (64 - s);
-    }
-  else if (s == 64)
-    {
-      i = p3l;
-      a = p2l;
-    }
-  else
-    { /* s > 64 */
-      i = p3l << (s - 64) | p2l >> (128 - s);
-      a = p2l << (s - 64) | p1l >> (128 - s);
-    }
-  int sgn = u;
-  sgn >>= 31;
-  int64_t sm = a >> 63;
-  i -= sm;
-  double z = (a ^ sgn) * 0x1p-64;
-  i = (i ^ sgn) - sgn;
-  *q = i;
-  return z;
-}
-
-static inline double
-rltl (float z, int *q)
-{
-  double x = z;
-  double idl = -0x1.b1bbead603d8bp-29 * x, idh = 0x1.45f306ep+2 * x,
-	 id = roundeven_finite (idh);
-  *q = asuint64 (0x1.8p52 + id);
-  return (idh - id) + idl;
-}
-
-static inline double
-rltl0 (double x, int *q)
-{
-  double idh = 0x1.45f306dc9c883p+2 * x, id = roundeven_finite (idh);
-  *q = asuint64 (0x1.8p52 + id);
-  return idh - id;
-}
-
-static inline float
-add_sign (float x, float rh, float rl)
-{
-  float sgn = copysignf (1.0f, x);
-  return sgn * rh + sgn * rl;
-}
-
 static float __attribute__ ((noinline))
 as_sinf_database (float x, double r)
 {
@@ -133,7 +66,7 @@ as_sinf_big (float x)
       return __math_invalidf (x);
     }
   int ia;
-  double z = rbig (t, &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]);
diff --git a/sysdeps/ieee754/flt-32/s_tanf.c b/sysdeps/ieee754/flt-32/s_tanf.c
index 5ee1d6f35e..d14c1f11bd 100644
--- a/sysdeps/ieee754/flt-32/s_tanf.c
+++ b/sysdeps/ieee754/flt-32/s_tanf.c
@@ -1,6 +1,6 @@
 /* Correctly-rounded tangent of binary32 value.
 
-Copyright (c) 2022-2024 Alexei Sibidanov.
+Copyright (c) 2022-2026 Alexei Sibidanov.
 
 The original version of this file was copied from the CORE-MATH
 project (file src/binary32/tan/tanf.c, revision 59d21d7).
@@ -27,13 +27,12 @@ SOFTWARE.
 #include <array_length.h>
 #include <stdint.h>
 #include <libm-alias-float.h>
-#include "math_config.h"
-#include <math_uint128.h>
+#include <s_sincosf_common.h>
 
 /* argument reduction
    for |z| < 2^28, return r such that 2/pi*x = q + r  */
-static inline double
-rltl (float z, int *q)
+static __always_inline double
+rltl_tanf (float z, int *q)
 {
   double x = z;
   double idl = -0x1.b1bbead603d8bp-32 * x;
@@ -43,58 +42,6 @@ rltl (float z, int *q)
   return (idh - id) + idl;
 }
 
-/* argument reduction
-   same as rltl, but for |x| >= 2^28  */
-static double __attribute__ ((noinline))
-rbig (uint32_t u, int *q)
-{
-  static const uint64_t ipi[] =
-    {
-      0xfe5163abdebbc562, 0xdb6295993c439041,
-      0xfc2757d1f534ddc0, 0xa2f9836e4e441529
-    };
-  int e = (u >> 23) & 0xff, i;
-  uint64_t m = (u & (~0u >> 9)) | 1 << 23;
-  u128 p0 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[0]));
-  u128 p1 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[1]));
-  p1 = u128_add (p1, u128_rshift (p0, 64));
-  u128 p2 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[2]));
-  p2 = u128_add (p2, u128_rshift (p1, 64));
-  u128 p3 = u128_mul (u128_from_u64 (m), u128_from_u64 (ipi[3]));
-  p3 = u128_add (p3, u128_rshift (p2, 64));
-  uint64_t p3h = u128_high (p3);
-  uint64_t p3l = u128_low (p3);
-  uint64_t p2l = u128_low (p2);
-  uint64_t p1l = u128_low (p1);
-  int64_t a;
-  int k = e - 127, s = k - 23;
-  /* in ctanf(), rbig() is called in the case 127+28 <= e < 0xff
-     thus 155 <= e <= 254, which yields 28 <= k <= 127 and 5 <= s <= 104 */
-  if (s < 64)
-    {
-      i = p3h << s | p3l >> (64 - s);
-      a = p3l << s | p2l >> (64 - s);
-    }
-  else if (s == 64)
-    {
-      i = p3l;
-      a = p2l;
-    }
-  else
-    { /* s > 64 */
-      i = p3l << (s - 64) | p2l >> (128 - s);
-      a = p2l << (s - 64) | p1l >> (128 - s);
-    }
-  int sgn = u;
-  sgn >>= 31;
-  int64_t sm = a >> 63;
-  i -= sm;
-  double z = (a ^ sgn) * 0x1p-64;
-  i = (i ^ sgn) - sgn;
-  *q = i;
-  return z;
-}
-
 float
 __tanf (float x)
 {
@@ -111,10 +58,10 @@ __tanf (float x)
 	  float x2 = x * x;
 	  return fmaf (x, 0x1.555556p-2f * x2, x);
 	}
-      z = rltl (x, &i);
+      z = rltl_tanf (x, &i);
     }
   else if (e < 0xff)
-    z = rbig (t, &i);
+    z = RBIG_TANF (t, &i);
   else
     {
       if (t << 9)
-- 
2.43.0