From patchwork Wed Mar 25 19:22:25 2026 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Adhemerval Zanella Netto X-Patchwork-Id: 132304 Return-Path: X-Original-To: patchwork@sourceware.org Delivered-To: patchwork@sourceware.org Received: from vm01.sourceware.org (localhost [127.0.0.1]) by sourceware.org (Postfix) with ESMTP id 267264B9DB55 for ; Wed, 25 Mar 2026 19:31:27 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 267264B9DB55 Authentication-Results: sourceware.org; dkim=pass (2048-bit key, unprotected) header.d=linaro.org header.i=@linaro.org header.a=rsa-sha256 header.s=google header.b=l0wIBn4G X-Original-To: libc-alpha@sourceware.org Delivered-To: libc-alpha@sourceware.org Received: from mail-vk1-xa36.google.com (mail-vk1-xa36.google.com [IPv6:2607:f8b0:4864:20::a36]) by sourceware.org (Postfix) with ESMTPS id C7BF24BB590E for ; Wed, 25 Mar 2026 19:24:22 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.2 sourceware.org C7BF24BB590E Authentication-Results: sourceware.org; dmarc=pass (p=none dis=none) header.from=linaro.org Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=linaro.org ARC-Filter: OpenARC Filter v1.0.0 sourceware.org C7BF24BB590E Authentication-Results: server2.sourceware.org; arc=none smtp.remote-ip=2607:f8b0:4864:20::a36 ARC-Seal: i=1; a=rsa-sha256; d=sourceware.org; s=key; t=1774466662; cv=none; b=M3/sI5sLRuTJIbaCCAETa+xvW600daq9KCho1/RptW+gUj0/S4hOcTZpixboFqNEQExGwhrlo2OEKF8iXV4t/GSm1SlPgPM0imVOj6wwU7qPp7LPg/G6MDIKNAloYaUxptVvgTW+3Tab3h76h3995/3dXM82pEbiljNUVs9xg7E= ARC-Message-Signature: i=1; a=rsa-sha256; d=sourceware.org; s=key; t=1774466662; c=relaxed/simple; bh=VoVpS9vd5ZHbsrxkIIi98zxb6cJeOPalQjtjKBNHI5U=; h=DKIM-Signature:From:To:Subject:Date:Message-ID:MIME-Version; b=a5f8TpCU5u9XCvOyAES3stuDxFQY2ZRvfX1VqzNPfiUsA4j5ugOClmu6s3dw1S1xWFo82MsrIpkf1joPBTLenpkd/KXSk6TXh44WBSQqmQsjTkrBqLyAj5Wpw187zve5qaoZNWbshhIaueKK63WylUO7PzZ7PcZOlzbLNnbsTKY= ARC-Authentication-Results: i=1; server2.sourceware.org DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org C7BF24BB590E Received: by mail-vk1-xa36.google.com with SMTP id 71dfb90a1353d-56ce1384618so200489e0c.1 for ; Wed, 25 Mar 2026 12:24:22 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=linaro.org; s=google; t=1774466662; x=1775071462; darn=sourceware.org; h=content-transfer-encoding:mime-version:references:in-reply-to :message-id:date:subject:cc:to:from:from:to:cc:subject:date :message-id:reply-to; bh=mbZwTzJo7Eo53aszNM2W6sjlbFFt4Lxr32/w2rjQeL4=; b=l0wIBn4GIjtOsXEKmjEVirnF+x454Eh+Bm9/sgSiI0OllUpE7Uhu/A+E9qkJ9vTkfG hXOQkry415rcMgchERN5sQIDgxGfhkjo973bsKt51wFx44+ADR5EEkh+XxfaqRrsww9g xenKMTseLOjeNyHlMWMzKm/CrTT65caq2V73XLaGMI79MIrH+P68ocVDNldDThHWOW4W pPVM6AyPew0TiHrPVIjJYDyXBZhUQmtbtgI6i3vVT18G3HoUNr0QjoMLOvZV0KefCSiQ lIhD/8PrPQAU1cKESZJsAgW1n9pf2y5KK/hjR5zpHS6PfKsHBVuiLtmA0iC2vcosZeV0 GiUw== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20251104; t=1774466662; x=1775071462; h=content-transfer-encoding:mime-version:references:in-reply-to :message-id:date:subject:cc:to:from:x-gm-gg:x-gm-message-state:from :to:cc:subject:date:message-id:reply-to; bh=mbZwTzJo7Eo53aszNM2W6sjlbFFt4Lxr32/w2rjQeL4=; b=Fh8tCze4Zs5cjpBM/BrSPJSHYxz4EvTsJAYwqPt7MGda4uX8HmHYKJ8kyAejKe1PLt JlPTYzFzaoikQMM9jLYMmYsQDM0NIGYf2orJfS7qnOWcTVhph31JCSiQkXC66LGCw7a6 CPZpMVKXjH9N+aQLnF+DChLGr804hTfFkX0coxSy+qVYIV/91q8GFBs4ql9Qsa/Pecgh i7saU9lo7jfPoEjNg3qnv0Z1Vu/wPNZA9lU/oz+cp0o/9diC46wmdFUtOVwGh6dRxny2 1pPhvV8ZtZ9QdR6FEk555QBE8agSBeaELUleS71F0bvnvXnXVho1FqESztoeiaPmlqd4 ZU1A== X-Gm-Message-State: AOJu0YxnlxniFQ2l9qjMP7vAKoxvcYA7b39S/erEs+c7NYsjo4xUJsUT zDzGV4Q2ngPgj2FbGhAQ6OXgwa50x1P1esC4JwjoCM6qdP5FDc3ivz6lcKu30PTRnnHKSQIjxxd rxfV4 X-Gm-Gg: ATEYQzx28IK+w1f/kBSfkW0ThcE2xTifmyRHjNmw9FJoFzrD8BwZv3uGH8M20j8qSH3 IrPLxLkR83rnBFmAxVb+OtCBc0lavNrnIec/0FU2yWKFXIZ6cveH1nOwj5/+evFQEXpCFNT2LwT p05O1VoLAjQM3F2+n9mfNG9CqzLG2Ng9MEXflEeM7C46IbgaIIaQu8SFN+v2HBkwl+33dp7pr0X j08gr1RuQmB2rY5v/aJl5FnJFOzVlRnM+bnulWimTjNgvQsorVr36JEdYut/u1NgbMLKieAB5hD LWsXxjuLS03wo1VCp7MBPASXDTT/ZO7UtHhTL3odjA7RyoNwE8iIc1hvut3/1fh3P0gDgdNzdWN A43uHgRl9ENmBj15Td/4FNHQ9PYjIqpvSh7ch70OhT7gaXvbINu5OvU0OmEy/WJxVzZ9KPGEUfp F4vaqXvdhF+83z3G2n3gc0wgPh32w6ni0+EEc= X-Received: by 2002:a05:6122:3d03:b0:563:7062:2a75 with SMTP id 71dfb90a1353d-56d21f9d5abmr2637692e0c.8.1774466661649; Wed, 25 Mar 2026 12:24:21 -0700 (PDT) Received: from mandiga.. ([2804:1b3:a7c1:90ea:f31d:ca7d:f8dd:c20b]) by smtp.gmail.com with ESMTPSA id 71dfb90a1353d-56d31d394fcsm1040021e0c.12.2026.03.25.12.24.19 (version=TLS1_3 cipher=TLS_AES_256_GCM_SHA384 bits=256/256); Wed, 25 Mar 2026 12:24:21 -0700 (PDT) From: Adhemerval Zanella To: libc-alpha@sourceware.org Cc: Paul Zimmermann , DJ Delorie Subject: [PATCH 7/8] math: Consolidate common definitions for cosf/sinf/tanf/sincosf Date: Wed, 25 Mar 2026 16:22:25 -0300 Message-ID: <20260325192357.1284741-8-adhemerval.zanella@linaro.org> X-Mailer: git-send-email 2.43.0 In-Reply-To: <20260325192357.1284741-1-adhemerval.zanella@linaro.org> References: <20260325192357.1284741-1-adhemerval.zanella@linaro.org> MIME-Version: 1.0 X-Spam-Status: No, score=-12.5 required=5.0 tests=BAYES_00, DKIM_SIGNED, DKIM_VALID, DKIM_VALID_AU, DKIM_VALID_EF, GIT_PATCH_0, RCVD_IN_DNSWL_NONE, SPF_HELO_NONE, SPF_PASS, TXREP, URIBL_BLOCKED autolearn=ham autolearn_force=no version=3.4.6 X-Spam-Checker-Version: SpamAssassin 3.4.6 (2021-04-09) on sourceware.org X-BeenThere: libc-alpha@sourceware.org X-Mailman-Version: 2.1.30 Precedence: list List-Id: Libc-alpha mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: libc-alpha-bounces~patchwork=sourceware.org@sourceware.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 #include "math_config.h" #include +#include #include -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 #include "math_config.h" #include +#include #include #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 +#include +#include + +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 +#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 #include "math_config.h" #include +#include #include #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 #include #include -#include "math_config.h" -#include +#include /* 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)