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 1106 invoked from network); 18 Aug 2023 21:17:53 -0000 Received: from second.openwall.net (193.110.157.125) by inbox.vuxu.org with ESMTPUTF8; 18 Aug 2023 21:17:53 -0000 Received: (qmail 17955 invoked by uid 550); 18 Aug 2023 21:17:50 -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 17915 invoked from network); 18 Aug 2023 21:17:49 -0000 Date: Fri, 18 Aug 2023 23:17:38 +0200 From: Szabolcs Nagy To: Paul Zimmermann Cc: musl@lists.openwall.com Message-ID: <20230818211738.GF3448312@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 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) powl used >= LDBL_MAX as infinity check, but LDBL_MAX is finite, so this can cause wrong results e.g. powl(LDBL_MAX, 0.5) returned inf or powl(2, LDBL_MAX) returned inf without raising overflow. huge y values (close to LDBL_MAX) could cause intermediate results to overflow (computing y * log2(x) with more than long double precision) and e.g. powl(0.5, 0x1p16380L) or powl(10, 0x1p16380L) returned nan. this is fixed by handling huge y early since that always overflows or underflows. reported by Paul Zimmermann against expl10 (which uses powl). --- src/math/powl.c | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/math/powl.c b/src/math/powl.c index 5b6da07b..6f64ea71 100644 --- a/src/math/powl.c +++ b/src/math/powl.c @@ -212,25 +212,33 @@ long double powl(long double x, long double y) } if (x == 1.0) return 1.0; /* 1**y = 1, even if y is nan */ - if (x == -1.0 && !isfinite(y)) - return 1.0; /* -1**inf = 1 */ if (y == 0.0) return 1.0; /* x**0 = 1, even if x is nan */ if (y == 1.0) return x; - if (y >= LDBL_MAX) { - if (x > 1.0 || x < -1.0) - return INFINITY; - if (x != 0.0) - return 0.0; - } - if (y <= -LDBL_MAX) { - if (x > 1.0 || x < -1.0) + /* if y*log2(x) < log2(LDBL_TRUE_MIN)-1 then x^y uflows to 0 + if y*log2(x) > -log2(LDBL_TRUE_MIN)+1 > LDBL_MAX_EXP then x^y oflows + if |x|!=1 then |log2(x)| > |log(x)| > LDBL_EPSILON/2 so + x^y oflows/uflows if |y|*LDBL_EPSILON/2 > -log2(LDBL_TRUE_MIN)+1 */ + if (fabsl(y) > 2*(-LDBL_MIN_EXP+LDBL_MANT_DIG+1)/LDBL_EPSILON) { + /* y is not an odd int */ + if (x == -1.0) + return 1.0; + if (y == INFINITY) { + if (x > 1.0 || x < -1.0) + return INFINITY; return 0.0; - if (x != 0.0 || y == -INFINITY) + } + if (y == -INFINITY) { + if (x > 1.0 || x < -1.0) + return 0.0; return INFINITY; + } + if ((x > 1.0 || x < -1.0) == (y > 0)) + return huge * huge; + return twom10000 * twom10000; } - if (x >= LDBL_MAX) { + if (x == INFINITY) { if (y > 0.0) return INFINITY; return 0.0; @@ -253,7 +261,7 @@ long double powl(long double x, long double y) yoddint = 1; } - if (x <= -LDBL_MAX) { + if (x == -INFINITY) { if (y > 0.0) { if (yoddint) return -INFINITY; -- 2.41.0