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.1 required=5.0 tests=DKIM_INVALID,DKIM_SIGNED, 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 20595 invoked from network); 19 Aug 2023 06:14:43 -0000 Received: from second.openwall.net (193.110.157.125) by inbox.vuxu.org with ESMTPUTF8; 19 Aug 2023 06:14:43 -0000 Received: (qmail 19491 invoked by uid 550); 19 Aug 2023 06:14:41 -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 19456 invoked from network); 19 Aug 2023 06:14:40 -0000 DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=inria.fr; s=dc; h=date:message-id:from:to:cc:in-reply-to:subject: references; bh=+6ia83pa1XHvFLdrJueXt2fCe2FytH/sf6d0tVC+Kzc=; b=cP3RY2j7K+Y0a+iytqMTY00FjnLKsEZQg7TqtM+KsF5v6ZBADjK3zlU3 5LmZJuAMZf45UJDf8t7skH6ufuDO3n+Jte+LYEZlubYF/ZeBA4cGvJn0d 1EMZNvZvM13N/eHDCWkL4FVPwbKeMolyS6Lu4Pn2z5MoncV2CzNdmVpEa c=; Authentication-Results: mail3-relais-sop.national.inria.fr; dkim=none (message not signed) header.i=none; spf=SoftFail smtp.mailfrom=Paul.Zimmermann@inria.fr; spf=None smtp.helo=postmaster@coriandre Received-SPF: SoftFail (mail3-relais-sop.national.inria.fr: domain of Paul.Zimmermann@inria.fr is inclined to not designate 152.81.9.227 as permitted sender) identity=mailfrom; client-ip=152.81.9.227; receiver=mail3-relais-sop.national.inria.fr; envelope-from="Paul.Zimmermann@inria.fr"; x-sender="Paul.Zimmermann@inria.fr"; x-conformance=spf_only; x-record-type="v=spf1"; x-record-text="v=spf1 ip4:128.93.142.0/24 ip4:192.134.164.0/24 ip4:128.93.162.160 ip4:89.107.174.7 mx ~all" Received-SPF: None (mail3-relais-sop.national.inria.fr: no sender authenticity information available from domain of postmaster@coriandre) identity=helo; client-ip=152.81.9.227; receiver=mail3-relais-sop.national.inria.fr; envelope-from="Paul.Zimmermann@inria.fr"; x-sender="postmaster@coriandre"; x-conformance=spf_only X-IronPort-AV: E=Sophos;i="6.01,185,1684792800"; d="scan'208";a="63764348" Date: Sat, 19 Aug 2023 08:14:28 +0200 Message-Id: From: Paul Zimmermann To: Szabolcs Nagy Cc: musl@lists.openwall.com In-Reply-To: <20230818211738.GF3448312@port70.net> (message from Szabolcs Nagy on Fri, 18 Aug 2023 23:17:38 +0200) References: <20230818211738.GF3448312@port70.net> Subject: [musl] Re: [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small) Dear Szabolcs, looks good to me, I now get a largest error of 40.1 ulps for exp10l: NEW exp10 0 0 0xd.4135bbaf96d9716p+8l [40] [40.1] 40.0525 40.05240557189884 Paul > Date: Fri, 18 Aug 2023 23:17:38 +0200 > From: Szabolcs Nagy > Cc: musl@lists.openwall.com > > 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 > >