From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from elaine.keithp.com (home.keithp.com [63.227.221.253]) by sourceware.org (Postfix) with ESMTPS id 055AC3858D35 for ; Sat, 1 Aug 2020 22:40:35 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 055AC3858D35 Received: from localhost (localhost [127.0.0.1]) by elaine.keithp.com (Postfix) with ESMTP id 30A7C3F2D17A; Sat, 1 Aug 2020 15:40:33 -0700 (PDT) X-Virus-Scanned: Debian amavisd-new at keithp.com Received: from elaine.keithp.com ([127.0.0.1]) by localhost (elaine.keithp.com [127.0.0.1]) (amavisd-new, port 10024) with LMTP id umQk1OyIbcBp; Sat, 1 Aug 2020 15:40:30 -0700 (PDT) Received: from keithp.com (koto.keithp.com [10.0.0.2]) by elaine.keithp.com (Postfix) with ESMTPSA id 6109C3F2D179; Sat, 1 Aug 2020 15:40:30 -0700 (PDT) Received: by keithp.com (Postfix, from userid 1000) id 47B8F15822C4; Sat, 1 Aug 2020 15:40:29 -0700 (PDT) From: "Keith Packard" To: Szabolcs Nagy Cc: newlib@sourceware.org Subject: Re: [PATCH] libm/math: Use __math_xflow in obsolete math code In-Reply-To: <20200731151941.GC5387@arm.com> References: <20200730234236.1282460-1-keithp@keithp.com> <20200731151941.GC5387@arm.com> Date: Sat, 01 Aug 2020 15:40:28 -0700 Message-ID: <877duihtwj.fsf@keithp.com> MIME-Version: 1.0 Content-Type: multipart/signed; boundary="==-=-="; micalg=pgp-sha256; protocol="application/pgp-signature" X-Spam-Status: No, score=-11.6 required=5.0 tests=BAYES_00, DKIM_SIGNED, DKIM_VALID, DKIM_VALID_AU, DKIM_VALID_EF, GIT_PATCH_0, 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: newlib@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Newlib mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 01 Aug 2020 22:40:38 -0000 --==-=-= Content-Type: multipart/mixed; boundary="=-=-=" --=-=-= Content-Type: text/plain Szabolcs Nagy writes: > this looks good. Thanks. > note1: i used c99 code when i wrote the > new math code and currently the old math > code can be compiled with older compilers, > this change might prevent that. I don't really have any way to test that as I'm using current compilers for this work. > note2: xflow may also try to set errno, > depending on the WANT_ERRNO setting (which > is on by default). I've got a separate patch ready which sets WANT_ERRNO based on _IEEE_LIBM so that both halves of the math library agree on whether errno should be set. That's sufficient for my uses (where I always set _IEEE_LIBM, and so never set errno) but not going to support the old distinction between the __ieee754 functions and the posix functions. I hadn't considered that issue. In my environment, I'm not enabling errno support. I can think of a couple of possible fixes: 1) Copy the __math_xflow functions to libm/math and remove errno handling so that the existing wrapper functions control how errno is managed. 2) Remove the __ieee754 entry points from the library so that the old library only offers the regular API, and the _IEEE_LIBM configuration changes those functions to not modify errno. Linking the new math code so that it also used the same configuration value would keep everyone in sync and offer matching errno behavior. Because the new math code uses double-precision arithmetic, I'm using the old math code on devices without double-precision hardware, like ARM processors that don't have __ARM_FP & 0x8. This includes pretty much all of the processors I build products out of. > in that case note that some compilers may model math functions as pure > (no sideeffect) and apply optimizations accrodingly, this may create > unexpected errno clobbers if math functions actually set errno. > (-fno-math-errno behaviour in gcc) but this affects other libm > implementations too and so far we haven't run into issues. It would be kinda cool if we could drive all of the errno support based on this compiler flag and eliminate the library configuration setting. We'd want this to take effect while building the application instead of at runtime. > i think e.g. return __math_uflowf(s == -1.0f) > would work better (but it's better to set a > sign flag where s is set). You're right -- here's an updated patch which uses (s < 0) as that doesn't require an extra constant. --=-=-= Content-Type: text/x-diff Content-Disposition: inline; filename=0001-libm-math-Use-__math_xflow-in-obsolete-math-code-v2.patch Content-Transfer-Encoding: quoted-printable From=20259c26b2a6697a813d5cd41923079eb873321b3b Mon Sep 17 00:00:00 2001 From: Keith Packard Date: Thu, 30 Jul 2020 16:41:05 -0700 Subject: [PATCH] libm/math: Use __math_xflow in obsolete math code [v2] C compilers may fold const values at compile time, so expressions which try to elicit underflow/overflow by performing simple arithemetic on suitable values will not generate the required exceptions. Work around this by replacing code which does these arithmetic operations with calls to the existing __math_xflow functions that are designed to do this correctly. Signed-off-by: Keith Packard =2D--- v2: libm/math: Pass sign to __math_xflow instead of muliplying result =2D-- newlib/libm/common/math_errf.c | 2 +- newlib/libm/math/e_cosh.c | 9 +++++---- newlib/libm/math/e_exp.c | 5 +++-- newlib/libm/math/e_pow.c | 18 ++++++++---------- newlib/libm/math/ef_cosh.c | 7 ++++--- newlib/libm/math/ef_exp.c | 5 +++-- newlib/libm/math/ef_pow.c | 14 ++++++-------- newlib/libm/math/s_erf.c | 3 ++- newlib/libm/math/sf_erf.c | 3 ++- 9 files changed, 34 insertions(+), 32 deletions(-) diff --git a/newlib/libm/common/math_errf.c b/newlib/libm/common/math_errf.c index 53c68b1cf..bb8273b8d 100644 =2D-- a/newlib/libm/common/math_errf.c +++ b/newlib/libm/common/math_errf.c @@ -51,13 +51,13 @@ xflowf (uint32_t sign, float y) return with_errnof (y, ERANGE); } =20 =2D#if !__OBSOLETE_MATH HIDDEN float __math_uflowf (uint32_t sign) { return xflowf (sign, 0x1p-95f); } =20 +#if !__OBSOLETE_MATH #if WANT_ERRNO_UFLOW /* Underflows to zero in some non-nearest rounding mode, setting errno is valid even if the result is non-zero, but in the subnormal range. */ diff --git a/newlib/libm/math/e_cosh.c b/newlib/libm/math/e_cosh.c index a6310bd07..7b258ffef 100644 =2D-- a/newlib/libm/math/e_cosh.c +++ b/newlib/libm/math/e_cosh.c @@ -25,7 +25,7 @@ * 2 * 22 <=3D x <=3D lnovft : cosh(x) :=3D exp(x)/2=20 * lnovft <=3D x <=3D ln2ovft: cosh(x) :=3D exp(x/2)/2 * exp(x/2) =2D * ln2ovft < x : cosh(x) :=3D huge*huge (overflow) + * ln2ovft < x : cosh(x) :=3D overflow * * Special cases: * cosh(x) is |x| if x is +INF, -INF, or NaN. @@ -33,13 +33,14 @@ */ =20 #include "fdlibm.h" +#include "math_config.h" =20 #ifndef _DOUBLE_IS_32BITS =20 #ifdef __STDC__ =2Dstatic const double one =3D 1.0, half=3D0.5, huge =3D 1.0e300; +static const double one =3D 1.0, half=3D0.5; #else =2Dstatic double one =3D 1.0, half=3D0.5, huge =3D 1.0e300; +static double one =3D 1.0, half=3D0.5; #endif =20 #ifdef __STDC__ @@ -87,7 +88,7 @@ static double one =3D 1.0, half=3D0.5, huge =3D 1.0e300; } =20 /* |x| > overflowthresold, cosh(x) overflow */ =2D return huge*huge; + return __math_oflow(0); } =20 #endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/math/e_exp.c b/newlib/libm/math/e_exp.c index d23b1162b..ec26c2099 100644 =2D-- a/newlib/libm/math/e_exp.c +++ b/newlib/libm/math/e_exp.c @@ -75,6 +75,7 @@ */ =20 #include "fdlibm.h" +#include "math_config.h" #if __OBSOLETE_MATH =20 #ifndef _DOUBLE_IS_32BITS @@ -126,8 +127,8 @@ P5 =3D 4.13813679705723846039e-08; /* 0x3E663769, 0x= 72BEA4D0 */ return x+x; /* NaN */ else return (xsb=3D=3D0)? x:0.0; /* exp(+-inf)=3D{inf,0} */ } =2D if(x > o_threshold) return huge*huge; /* overflow */ =2D if(x < u_threshold) return twom1000*twom1000; /* underflow */ + if(x > o_threshold) return __math_oflow(0); /* overflow */ + if(x < u_threshold) return __math_uflow(0); /* underflow */ } =20 /* argument reduction */ diff --git a/newlib/libm/math/e_pow.c b/newlib/libm/math/e_pow.c index 4c450ec05..258cca8dd 100644 =2D-- a/newlib/libm/math/e_pow.c +++ b/newlib/libm/math/e_pow.c @@ -76,8 +76,6 @@ zero =3D 0.0, one =3D 1.0, two =3D 2.0, two53 =3D 9007199254740992.0, /* 0x43400000, 0x00000000 */ =2Dhuge =3D 1.0e300, =2Dtiny =3D 1.0e-300, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ L1 =3D 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ L2 =3D 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ @@ -197,12 +195,12 @@ ivln2_l =3D 1.92596299112661746887e-08; /* 0x3E54AE= 0B, 0xF85DDF44 =3D1/ln2 tail*/ /* |y| is huge */ if(iy>0x41e00000) { /* if |y| > 2**31 */ if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ =2D if(ix<=3D0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; =2D if(ix>=3D0x3ff00000) return (hy>0)? huge*huge:tiny*tiny; + if(ix<=3D0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0); + if(ix>=3D0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0); } /* over/underflow if x is not close to one */ =2D if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; =2D if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny; + if(ix<0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0); + if(ix>0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0); /* now |1-x| is tiny <=3D 2**-20, suffice to compute=20 log(x) by x-x^2/2+x^3/3-x^4/4 */ t =3D ax-1; /* t has 20 trailing zeros */ @@ -275,15 +273,15 @@ ivln2_l =3D 1.92596299112661746887e-08; /* 0x3E54AE= 0B, 0xF85DDF44 =3D1/ln2 tail*/ EXTRACT_WORDS(j,i,z); if (j>=3D0x40900000) { /* z >=3D 1024 */ if(((j-0x40900000)|i)!=3D0) /* if z > 1024 */ =2D return s*huge*huge; /* overflow */ + return __math_oflow(s<0); /* overflow */ else { =2D if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ + if(p_l+ovt>z-p_h) return __math_oflow(s<0); /* overflow */ } } else if((j&0x7fffffff)>=3D0x4090cc00 ) { /* z <=3D -1075 */ if(((j-0xc090cc00)|i)!=3D0) /* z < -1075 */ =2D return s*tiny*tiny; /* underflow */ + return __math_uflow(s<0); /* underflow */ else { =2D if(p_l<=3Dz-p_h) return s*tiny*tiny; /* underflow */ + if(p_l<=3Dz-p_h) return __math_uflow(s<0); /* underflow */ } } /* diff --git a/newlib/libm/math/ef_cosh.c b/newlib/libm/math/ef_cosh.c index bdce61a00..5690bd7a4 100644 =2D-- a/newlib/libm/math/ef_cosh.c +++ b/newlib/libm/math/ef_cosh.c @@ -14,15 +14,16 @@ */ =20 #include "fdlibm.h" +#include "math_config.h" =20 #ifdef __v810__ #define const #endif =20 #ifdef __STDC__ =2Dstatic const float one =3D 1.0, half=3D0.5, huge =3D 1.0e30; +static const float one =3D 1.0, half=3D0.5; #else =2Dstatic float one =3D 1.0, half=3D0.5, huge =3D 1.0e30; +static float one =3D 1.0, half=3D0.5; #endif =20 #ifdef __STDC__ @@ -67,5 +68,5 @@ static float one =3D 1.0, half=3D0.5, huge =3D 1.0e30; } =20 /* |x| > overflowthresold, cosh(x) overflow */ =2D return huge*huge; + return __math_oflowf(0); } diff --git a/newlib/libm/math/ef_exp.c b/newlib/libm/math/ef_exp.c index fb3e2ffe6..0cd0c00b3 100644 =2D-- a/newlib/libm/math/ef_exp.c +++ b/newlib/libm/math/ef_exp.c @@ -14,6 +14,7 @@ */ =20 #include "fdlibm.h" +#include "math_config.h" =20 #if __OBSOLETE_MATH #ifdef __v810__ @@ -61,9 +62,9 @@ P5 =3D 4.1381369442e-08; /* 0x3331bb4c */ if(FLT_UWORD_IS_INFINITE(hx)) return (xsb=3D=3D0)? x:0.0; /* exp(+-inf)=3D{inf,0} */ if(sx > FLT_UWORD_LOG_MAX) =2D return huge*huge; /* overflow */ + return __math_oflowf(0); /* overflow */ if(sx < 0 && hx > FLT_UWORD_LOG_MIN) =2D return twom100*twom100; /* underflow */ + return __math_uflow(0); /* underflow */ =20=09 /* argument reduction */ if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */=20 diff --git a/newlib/libm/math/ef_pow.c b/newlib/libm/math/ef_pow.c index d4ea4c5e8..e3579f071 100644 =2D-- a/newlib/libm/math/ef_pow.c +++ b/newlib/libm/math/ef_pow.c @@ -33,8 +33,6 @@ zero =3D 0.0, one =3D 1.0, two =3D 2.0, two24 =3D 16777216.0, /* 0x4b800000 */ =2Dhuge =3D 1.0e30, =2Dtiny =3D 1.0e-30, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ L1 =3D 6.0000002384e-01, /* 0x3f19999a */ L2 =3D 4.2857143283e-01, /* 0x3edb6db7 */ @@ -140,8 +138,8 @@ ivln2_l =3D 7.0526075433e-06; /* 0x36eca570 =3D1/ln2 = tail*/ /* |y| is huge */ if(iy>0x4d000000) { /* if |y| > 2**27 */ /* over/underflow if x is not close to one */ =2D if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny; =2D if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny; + if(ix<0x3f7ffff8) return (hy<0)? __math_oflowf(0):__math_uflowf(0); + if(ix>0x3f800007) return (hy>0)? __math_oflowf(0):__math_uflowf(0); /* now |1-x| is tiny <=3D 2**-20, suffice to compute=20 log(x) by x-x^2/2+x^3/3-x^4/4 */ t =3D ax-1; /* t has 20 trailing zeros */ @@ -219,14 +217,14 @@ ivln2_l =3D 7.0526075433e-06; /* 0x36eca570 =3D1/ln= 2 tail*/ i =3D j&0x7fffffff; if (j>0) { if (i>FLT_UWORD_EXP_MAX) =2D return s*huge*huge; /* overflow */ + return __math_oflowf(s<0); /* overflow */ else if (i=3D=3DFLT_UWORD_EXP_MAX) =2D if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ + if(p_l+ovt>z-p_h) return __math_oflowf(s<0); /* overflow */ } else { if (i>FLT_UWORD_EXP_MIN) =2D return s*tiny*tiny; /* underflow */ + return __math_uflowf(s<0); /* underflow */ else if (i=3D=3DFLT_UWORD_EXP_MIN) =2D if(p_l<=3Dz-p_h) return s*tiny*tiny; /* underflow */ + if(p_l<=3Dz-p_h) return __math_uflowf(s<0); /* underflow */ } /* * compute 2**(p_h+p_l) diff --git a/newlib/libm/math/s_erf.c b/newlib/libm/math/s_erf.c index eb288fc73..9e2333c11 100644 =2D-- a/newlib/libm/math/s_erf.c +++ b/newlib/libm/math/s_erf.c @@ -152,6 +152,7 @@ PORTABILITY =20 =20 #include "fdlibm.h" +#include "math_config.h" =20 #ifndef _DOUBLE_IS_32BITS =20 @@ -352,7 +353,7 @@ sb7 =3D -2.24409524465858183362e+01; /* 0xC03670E2, 0x= 42712D62 */ __ieee754_exp((z-x)*(z+x)+R/S); if(hx>0) return r/x; else return two-r/x; } else { =2D if(hx>0) return tiny*tiny; else return two-tiny; + if(hx>0) return __math_uflow(0); else return two-tiny; } } =20 diff --git a/newlib/libm/math/sf_erf.c b/newlib/libm/math/sf_erf.c index 0329c60fa..f3d0de97a 100644 =2D-- a/newlib/libm/math/sf_erf.c +++ b/newlib/libm/math/sf_erf.c @@ -14,6 +14,7 @@ */ =20 #include "fdlibm.h" +#include "math_config.h" =20 #ifdef __v810__ #define const @@ -217,7 +218,7 @@ sb7 =3D -2.2440952301e+01; /* 0xc1b38712 */ __ieee754_expf((z-x)*(z+x)+R/S); if(hx>0) return r/x; else return two-r/x; } else { =2D if(hx>0) return tiny*tiny; else return two-tiny; + if(hx>0) return __math_uflow(0); else return two-tiny; } } =20 =2D-=20 2.28.0 --=-=-= Content-Type: text/plain Content-Transfer-Encoding: quoted-printable =2D-=20 =2Dkeith --=-=-=-- --==-=-= Content-Type: application/pgp-signature; name="signature.asc" -----BEGIN PGP SIGNATURE----- iQIzBAEBCAAdFiEEw4O3eCVWE9/bQJ2R2yIaaQAAABEFAl8l71wACgkQ2yIaaQAA ABF2JQ//T5ktDRGVR6JKWxl4COhbzXPBZE7de4FoCYLjDDvis0mcJymCYtxfPhWx 0mQ4YTGOGljx87AotFJYjw1J0GJ0I5NjsReTmekvFlpc4ST9GhlDLe+IgSDBAx4x VmbZ685tEPaYR1CqfMRvmg2EuoCeLsTpApr6CUpVKevcwH1U3u9QEoDtlAY4GKtD WZkVQB9vt3/EnVGzWg4OooV4pjJKve1YdqL0uwziSxSGh6IWyTqsEnYq1vTjuuZS Hh5J/7uE6AFgkfjdk1g4jTPIrxa5LSniXn/Yv2ZZnjBklktInENhG/rzXpsADitz yEqFOyRh1AY8U4C67+kzGo1Ud6hhGb0eGY3KQ+69hzt9jyigt7bjiPiG3/hUBwFU AmYk/k8kLHn7NA2gxMmVgl64+24P70C024gRCEpc36EvIs0NLfQJFeOj6ThGJCGE uK7b7Q2JZgG5qltHrmzfHH6aluqao0uQuyi1j+Y6OtfOD4SKMVwzOhUOuuUdpdVz epHcn95/DdrnG+WDB4jHAFsnEon8FY9cGnJJh7fVnN8t9dMdI3Ta90Wwd/XQVfqC s+Nx1L96lydBNjJaPTsi+jdgKh1eoatozIWHCI4S/5QlskC1gShiFZoAWd7INuQQ bQwoL9srU4bjTEikpL7OpjN5bp4p1GQh9j/99Y34Ao8TZciP8oM= =uk5J -----END PGP SIGNATURE----- --==-=-=--