From patchwork Wed Mar 31 19:27:53 2021 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Paul Zimmermann X-Patchwork-Id: 42824 Return-Path: X-Original-To: patchwork@sourceware.org Delivered-To: patchwork@sourceware.org Received: from server2.sourceware.org (localhost [IPv6:::1]) by sourceware.org (Postfix) with ESMTP id 583F03857C74; Wed, 31 Mar 2021 19:28:24 +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 921153857C6F for ; Wed, 31 Mar 2021 19:28:21 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 921153857C6F 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:b8ThyKwGuNkcbqO0NNa+KrPw1b1zdoIgy1knxilNYDZeG/b1q+mFmvMH2RjozAsLUHY7ltyafIWGS3XQ9Zl6iLNhX4uKcQH6tAKTQr1KwpDlx1TbcRHW0s54+eNef7NlCNv2ZGIK7vrSxAWjCd4vzJ2m/cmT5Nv29HtmQQF0Z6wI1W4QYTqzKEFwSQVcbKBJcaa03NZNpDarZB0sAfiTO39tZYX+juyOsJrnZBIcbiRG1DWz X-IronPort-AV: E=Sophos;i="5.81,293,1610406000"; d="scan'208";a="500978380" Received: from tomate.loria.fr ([152.81.10.51]) by mail2-relais-roc.national.inria.fr with ESMTP/TLS/DHE-RSA-AES256-GCM-SHA384; 31 Mar 2021 21:27:55 +0200 Received: from zimmerma by tomate.loria.fr with local (Exim 4.94) (envelope-from ) id 1lRgVO-00BAJH-OV; Wed, 31 Mar 2021 21:27:54 +0200 From: Paul Zimmermann To: libc-alpha@sourceware.org Subject: [PATCH 2/2] added auxiliary routine reduce_aux Date: Wed, 31 Mar 2021 21:27:53 +0200 Message-Id: <20210331192753.2661194-1-Paul.Zimmermann@inria.fr> X-Mailer: git-send-email 2.30.2 MIME-Version: 1.0 X-Spam-Status: No, score=-10.1 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.h | 64 +++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 sysdeps/ieee754/flt-32/reduce_aux.h diff --git a/sysdeps/ieee754/flt-32/reduce_aux.h b/sysdeps/ieee754/flt-32/reduce_aux.h new file mode 100644 index 0000000000..394721a9fd --- /dev/null +++ b/sysdeps/ieee754/flt-32/reduce_aux.h @@ -0,0 +1,64 @@ +/* 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 + . */ + +#ifndef _MATH_REDUCE_AUX_H +#define _MATH_REDUCE_AUX_H + +#include +#include +#include + +/* 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; +} + +#endif