From patchwork Tue Feb 9 05:28:15 2021 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Paul Zimmermann X-Patchwork-Id: 1438114 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Authentication-Results: ozlabs.org; spf=pass (sender SPF authorized) smtp.mailfrom=sourceware.org (client-ip=8.43.85.97; helo=sourceware.org; envelope-from=libc-alpha-bounces@sourceware.org; receiver=) Received: from sourceware.org (unknown [8.43.85.97]) (using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) key-exchange X25519 server-signature RSA-PSS (4096 bits) server-digest SHA256) (No client certificate requested) by ozlabs.org (Postfix) with ESMTPS id 4DZWdS5tFcz9rx8 for ; Tue, 9 Feb 2021 16:28:28 +1100 (AEDT) Received: from server2.sourceware.org (localhost [IPv6:::1]) by sourceware.org (Postfix) with ESMTP id BB17A3861000; Tue, 9 Feb 2021 05:28:25 +0000 (GMT) X-Original-To: libc-alpha@sourceware.org Delivered-To: libc-alpha@sourceware.org Received: from mail2-relais-roc.national.inria.fr (mail2-relais-roc.national.inria.fr [192.134.164.83]) by sourceware.org (Postfix) with ESMTPS id 1D0D03861000 for ; Tue, 9 Feb 2021 05:28:22 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 1D0D03861000 Authentication-Results: sourceware.org; dmarc=none (p=none dis=none) header.from=inria.fr Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=Paul.Zimmermann@loria.fr X-IronPort-AV: E=Sophos;i="5.81,164,1610406000"; d="scan'208";a="491909311" Received: from tomate.loria.fr ([152.81.10.51]) by mail2-relais-roc.national.inria.fr with ESMTP/TLS/DHE-RSA-AES256-GCM-SHA384; 09 Feb 2021 06:28:21 +0100 Received: from zimmerma by tomate.loria.fr with local (Exim 4.94) (envelope-from ) id 1l9LZU-00GkJq-Cw; Tue, 09 Feb 2021 06:28:20 +0100 From: Paul Zimmermann To: libc-alpha@sourceware.org Subject: [PATCH 1/4] Add a routine for accurate reduction mod (2pi), for j0f/j1f/y0f/y1f. Date: Tue, 9 Feb 2021 06:28:15 +0100 Message-Id: <20210209052818.3991249-1-Paul.Zimmermann@inria.fr> X-Mailer: git-send-email 2.30.0 MIME-Version: 1.0 X-Spam-Status: No, score=-9.7 required=5.0 tests=BAYES_00, GIT_PATCH_0, HEADER_FROM_DIFFERENT_DOMAINS, KAM_DMARC_STATUS, KAM_SHORT, RCVD_IN_MSPIKE_H2, SPF_HELO_NONE, SPF_PASS, TXREP autolearn=ham autolearn_force=no version=3.4.2 X-Spam-Checker-Version: SpamAssassin 3.4.2 (2018-09-13) on server2.sourceware.org X-BeenThere: libc-alpha@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Libc-alpha mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: libc-alpha-bounces@sourceware.org Sender: "Libc-alpha" --- math/Makefile | 3 +- sysdeps/ieee754/flt-32/e_reduce_pi.c | 186 +++++++++++++++++++++++++++ sysdeps/ieee754/flt-32/e_reduce_pi.h | 19 +++ 3 files changed, 207 insertions(+), 1 deletion(-) create mode 100644 sysdeps/ieee754/flt-32/e_reduce_pi.c create mode 100644 sysdeps/ieee754/flt-32/e_reduce_pi.h diff --git a/math/Makefile b/math/Makefile index 687aa5d510..6f5dfcfd42 100644 --- a/math/Makefile +++ b/math/Makefile @@ -134,7 +134,8 @@ type-double-routines := branred doasin dosincos mpa mpatan2 \ # float support type-float-suffix := f type-float-routines := math_errf e_exp2f_data e_logf_data \ - e_log2f_data e_powf_log2_data s_sincosf_data + e_log2f_data e_powf_log2_data s_sincosf_data \ + e_reduce_pi # _Float128 support type-float128-suffix := f128 diff --git a/sysdeps/ieee754/flt-32/e_reduce_pi.c b/sysdeps/ieee754/flt-32/e_reduce_pi.c new file mode 100644 index 0000000000..2028d44817 --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_reduce_pi.c @@ -0,0 +1,186 @@ +/* Argument reduction modulo 2pi. + Copyright (C) 2021 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 + . */ + +#include "math_config.h" + +/* Argument reduction mod 2pi: T[i] ~ 2^i mod (2*pi). */ +static const double T[128] = { + 0x1p+0, + 0x2p+0, + -0x2.487ed5110b462p+0, + 0x1.b7812aeef4b9fp+0, + -0x2.d97c7f3321d24p+0, + 0x9.585d6aac7a1a8p-4, + 0x1.2b0bad558f435p+0, + 0x2.56175aab1e86ap+0, + -0x1.9c501fbace38dp+0, + 0x3.0fde959b6ed46p+0, + -0x2.8c1a9da2d9d3cp-4, + -0x5.18353b45b3a78p-4, + -0xa.306a768b674fp-4, + -0x1.460d4ed16ce9ep+0, + -0x2.8c1a9da2d9d3cp+0, + 0x1.304999cb579e8p+0, + 0x2.60933396af3dp+0, + -0x1.87586de3accc3p+0, + -0x3.0eb0dbc759986p+0, + 0x2.b1d1d8258155p-4, + 0x5.63a3b04b02aap-4, + 0xa.c74760960554p-4, + 0x1.58e8ec12c0aa8p+0, + 0x2.b1d1d8258155p+0, + -0xe.4db24c6089bf8p-4, + -0x1.c9b6498c1137fp+0, + 0x2.b51241f8e8d64p+0, + -0xd.e5a511f3999bp-4, + -0x1.bcb4a23e73336p+0, + 0x2.cf15909424df4p+0, + -0xa.a53b3e8c1877p-4, + -0x1.54a767d1830eep+0, + -0x2.a94ecfa3061dcp+0, + 0xf.5e135caff0a78p-4, + 0x1.ebc26b95fe14fp+0, + -0x2.70f9fde50f1c4p+0, + 0x1.668ad946ed0dap+0, + 0x2.cd15b28dda1b4p+0, + -0xa.e536ff5570fbp-4, + -0x1.5ca6dfeaae1f6p+0, + -0x2.b94dbfd55c3ecp+0, + 0xd.5e3556652c89p-4, + 0x1.abc6aacca5912p+0, + -0x2.f0f17f77c023ep+0, + 0x6.69bd6218afe64p-4, + 0xc.d37ac4315fcc8p-4, + 0x1.9a6f58862bf99p+0, + -0x3.13a02404b352ep+0, + 0x2.13e8d07a4a046p-4, + 0x4.27d1a0f49408cp-4, + 0x8.4fa341e928118p-4, + 0x1.09f4683d25023p+0, + 0x2.13e8d07a4a046p+0, + -0x2.20ad341c773d4p+0, + 0x2.07246cd81ccb8p+0, + -0x2.3a35fb60d1af2p+0, + 0x1.d412de4f67e7ep+0, + -0x2.a05918723b764p+0, + 0x1.07cca42c94599p+0, + 0x2.0f99485928b32p+0, + -0x2.294c445eb9dfep+0, + 0x1.f5e64c5397865p+0, + -0x2.5cb23c69dc396p+0, + 0x1.8f1a5c3d52d34p+0, + 0x3.1e34b87aa5a68p+0, + -0xc.15641bbff90bp-8, + -0x1.82ac8377ff216p-4, + -0x3.055906effe42cp-4, + -0x6.0ab20ddffc858p-4, + -0xc.15641bbff90bp-4, + -0x1.82ac8377ff216p+0, + -0x3.055906effe42cp+0, + 0x3.dccc7310ec09ap-4, + 0x7.b998e621d8134p-4, + 0xf.7331cc43b0268p-4, + 0x1.ee6639887604dp+0, + -0x2.6bb262001f3c6p+0, + 0x1.711a1110cccd4p+0, + 0x2.e2342221999a8p+0, + -0x8.41690cdd81128p-4, + -0x1.082d219bb0225p+0, + -0x2.105a43376044ap+0, + 0x2.27ca4ea24abcep+0, + -0x1.f8ea37cc75cc6p+0, + 0x2.56aa65781fad6p+0, + -0x1.9b2a0a20cbeb6p+0, + 0x3.122ac0cf736f6p+0, + -0x2.4295372246764p-4, + -0x4.852a6e448cec8p-4, + -0x9.0a54dc8919d9p-4, + -0x1.214a9b91233b2p+0, + -0x2.4295372246764p+0, + 0x1.c35466cc7e59ap+0, + -0x2.c1d607780e92ep+0, + 0xc.4d2c620ee206p-4, + 0x1.89a58c41dc40cp+0, + 0x3.134b1883b8818p+0, + -0x2.1e8a4099a4324p-4, + -0x4.3d14813348648p-4, + -0x8.7a29026690c9p-4, + -0x1.0f45204cd2192p+0, + -0x2.1e8a4099a4324p+0, + 0x2.0b6a53ddc2e18p+0, + -0x2.31aa2d558583p+0, + 0x1.e52a7a6600402p+0, + -0x2.7e29e0450ac5cp+0, + 0x1.4c2b1486f5ba8p+0, + 0x2.9856290deb75p+0, + -0x1.17d282f5345bfp+0, + -0x2.2fa505ea68b7ep+0, + 0x1.e934c93c39d64p+0, + -0x2.761542989799ap+0, + 0x1.5c544fdfdc12fp+0, + 0x2.b8a89fbfb825ep+0, + -0xd.72d95919afa58p-4, + -0x1.ae5b2b2335f4bp+0, + 0x2.ebc87eca9f5cap+0, + -0x7.0edd77bcc8ccp-4, + -0xe.1dbaef799198p-4, + -0x1.c3b75def3233p+0, + 0x2.c1101932a6e02p+0, + -0xc.65ea2abbd85fp-4, + -0x1.8cbd45577b0bep+0, + -0x3.197a8aaef617cp+0, + 0x1.589bfb31f1687p-4, + 0x2.b137f663e2d0ep-4, + 0x5.626fecc7c5a1cp-4, + 0xa.c4dfd98f8b438p-4 +}; + +/* Return h such that x - pi/4 mod (2*pi) ~ h, + with maximal error about 2^-45. */ +attribute_hidden double +__math_reduce_mod_twopi (double x) +{ + double sign = 1, h; + if (x < 0) + { + x = -x; + sign = -1; + } + h = 0; /* Invariant is x+h mod (2*pi). */ + while (x >= 4) + { + int i = ilogb (x); + x -= ldexp (1.0f, i); + h += T[i]; + } + /* Add the remainder x. */ + h += x; + /* Subtract pi/4. */ + double piover4 = 0xc.90fdaa22168cp-4; + h -= piover4; + /* Reduce mod 2*pi. */ + double twopi = 0x6.487ed5110b46p+0; + while (fabs (h) > twopi * 0.5) + { + if (h > 0) + h -= twopi; + else + h += twopi; + } + return sign * h; +} diff --git a/sysdeps/ieee754/flt-32/e_reduce_pi.h b/sysdeps/ieee754/flt-32/e_reduce_pi.h new file mode 100644 index 0000000000..93d6ce8764 --- /dev/null +++ b/sysdeps/ieee754/flt-32/e_reduce_pi.h @@ -0,0 +1,19 @@ +/* Argument reduction modulo 2pi. + Copyright (C) 2021 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 + . */ + +attribute_hidden double __math_reduce_mod_twopi (double);