From patchwork Fri Mar 19 15:06:18 2021 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Paul Zimmermann X-Patchwork-Id: 1455918 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 4F26g6211rz9sVS for ; Sat, 20 Mar 2021 02:06:42 +1100 (AEDT) Received: from server2.sourceware.org (localhost [IPv6:::1]) by sourceware.org (Postfix) with ESMTP id 24F00384A02E; Fri, 19 Mar 2021 15:06:39 +0000 (GMT) X-Original-To: libc-alpha@sourceware.org Delivered-To: libc-alpha@sourceware.org Received: from mail3-relais-sop.national.inria.fr (mail3-relais-sop.national.inria.fr [192.134.164.104]) by sourceware.org (Postfix) with ESMTPS id 1E4093857C61 for ; Fri, 19 Mar 2021 15:06:36 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 1E4093857C61 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 IronPort-HdrOrdr: A9a23:3k3e/qkWQfdMMj7s+CN6QsjeGI7pDfLV3DAbvn1ZSRFFG/GwvcrGppgm/DXzjyscX2xlpMuJP7OOTWiZ2Zl+54QQOrnKZniFhEKEJJxvhLGN/xTOACv7n9Q26Y5FU4xTTOL9FkJ7i8GS2njaL/8FzMOc+K6lwcfSpk0NcShQZ6tt7xh0B2+geyUceCB9CZ01GIH03LsjmxObZX8VYs6nb0NrY8H/obTw+a7OUFojDx4j5BLmt1OV1II= X-IronPort-AV: E=Sophos;i="5.81,262,1610406000"; d="scan'208";a="376263676" Received: from tomate.loria.fr ([152.81.10.51]) by mail3-relais-sop.national.inria.fr with ESMTP/TLS/DHE-RSA-AES256-GCM-SHA384; 19 Mar 2021 16:06:31 +0100 Received: from zimmerma by tomate.loria.fr with local (Exim 4.94) (envelope-from ) id 1lNGhq-0053z8-VE; Fri, 19 Mar 2021 16:06:31 +0100 From: Paul Zimmermann To: libc-alpha@sourceware.org Subject: [PATCH 1/9] Auxiliary function for reduction modulo 2*pi. Date: Fri, 19 Mar 2021 16:06:18 +0100 Message-Id: <20210319150626.1206905-1-Paul.Zimmermann@inria.fr> X-Mailer: git-send-email 2.30.2 MIME-Version: 1.0 X-Spam-Status: No, score=-10.0 required=5.0 tests=BAYES_00, GIT_PATCH_0, HEADER_FROM_DIFFERENT_DOMAINS, KAM_DMARC_STATUS, KAM_SHORT, RCVD_IN_MSPIKE_H3, RCVD_IN_MSPIKE_WL, 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" --- sysdeps/ieee754/flt-32/reduce_aux.c | 55 +++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 sysdeps/ieee754/flt-32/reduce_aux.c Reviewed-by: Adhemerval Zanella diff --git a/sysdeps/ieee754/flt-32/reduce_aux.c b/sysdeps/ieee754/flt-32/reduce_aux.c new file mode 100644 index 0000000000..412b4d22cb --- /dev/null +++ b/sysdeps/ieee754/flt-32/reduce_aux.c @@ -0,0 +1,55 @@ +/* Auxiliary routine for the Bessel functions (j0f, y0f, j1f, y1f). + 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 + . */ + +/* Return h and update n such that: + Now x - pi/4 - alpha = h + n*pi/2 mod (2*pi). */ +static inline double +reduce_aux (float x, int *n, double alpha) +{ + double h; + h = reduce_large (asuint (x), n); + /* Now |x| = h+n*pi/2 mod 2*pi. */ + /* Recover sign. */ + if (x < 0) + { + h = -h; + *n = -*n; + } + /* Subtract pi/4. */ + double piover2 = 0xc.90fdaa22168cp-3; + if (h >= 0) + h -= piover2 / 2; + else + { + h += piover2 / 2; + (*n) --; + } + /* Subtract alpha and reduce if needed mod pi/2. */ + h -= alpha; + if (h > piover2) + { + h -= piover2; + (*n) ++; + } + else if (h < -piover2) + { + h += piover2; + (*n) --; + } + return h; +}