From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on inbox.vuxu.org X-Spam-Level: X-Spam-Status: No, score=-3.3 required=5.0 tests=MAILING_LIST_MULTI, RCVD_IN_DNSWL_MED,RCVD_IN_MSPIKE_H3,RCVD_IN_MSPIKE_WL autolearn=ham autolearn_force=no version=3.4.4 Received: (qmail 1048 invoked from network); 18 Aug 2023 21:16:17 -0000 Received: from second.openwall.net (193.110.157.125) by inbox.vuxu.org with ESMTPUTF8; 18 Aug 2023 21:16:17 -0000 Received: (qmail 15629 invoked by uid 550); 18 Aug 2023 21:16:13 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Reply-To: musl@lists.openwall.com Received: (qmail 15594 invoked from network); 18 Aug 2023 21:16:12 -0000 Date: Fri, 18 Aug 2023 23:16:00 +0200 From: Szabolcs Nagy To: Paul Zimmermann Cc: musl@lists.openwall.com Message-ID: <20230818211600.GE3448312@port70.net> Mail-Followup-To: Paul Zimmermann , musl@lists.openwall.com References: MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: Subject: [musl] [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0 acosh(x) is nan for x < 1, but x < 0 cases were not handled specially and acoshl gave wrong result for some -0x1p32 < x < -2 values, e.g.: acoshl(-0x1p20) returned -inf, acoshl(-0x1.4p20) returned -0x1.db365758403aa9acp+0L, fixed by checking the sign bit and handling it specially. reported by Paul Zimmermann. --- src/math/acoshl.c | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/math/acoshl.c b/src/math/acoshl.c index 8d4b43f6..943cec17 100644 --- a/src/math/acoshl.c +++ b/src/math/acoshl.c @@ -10,14 +10,18 @@ long double acoshl(long double x) long double acoshl(long double x) { union ldshape u = {x}; - int e = u.i.se & 0x7fff; + int e = u.i.se; if (e < 0x3fff + 1) - /* |x| < 2, invalid if x < 1 or nan */ + /* 0 <= x < 2, invalid if x < 1 */ return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1))); if (e < 0x3fff + 32) - /* |x| < 0x1p32 */ + /* 2 <= x < 0x1p32 */ return logl(2*x - 1/(x+sqrtl(x*x-1))); + if (e & 0x8000) + /* x < 0 or x = -0, invalid */ + return (x - x) / (x - x); + /* 0x1p32 <= x or nan */ return logl(x) + 0.693147180559945309417232121458176568L; } #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 -- 2.41.0