@@ -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
@@ -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]);
@@ -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]);
new file mode 100644
@@ -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);
+}
new file mode 100644
@@ -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
@@ -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]);
@@ -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)