From mboxrd@z Thu Jan 1 00:00:00 1970 X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/6135 Path: news.gmane.org!not-for-mail From: Sergey Dmitrouk Newsgroups: gmane.linux.lib.musl.general Subject: [PATCH] Make musl math depend less on libgcc builtins Date: Thu, 11 Sep 2014 10:35:32 +0300 Message-ID: <20140911073532.GA3179@zx-spectrum> Reply-To: musl@lists.openwall.com NNTP-Posting-Host: plane.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="LZvS9be/3tNcYl/X" X-Trace: ger.gmane.org 1410420958 15901 80.91.229.3 (11 Sep 2014 07:35:58 GMT) X-Complaints-To: usenet@ger.gmane.org NNTP-Posting-Date: Thu, 11 Sep 2014 07:35:58 +0000 (UTC) To: Original-X-From: musl-return-6148-gllmg-musl=m.gmane.org@lists.openwall.com Thu Sep 11 09:35:52 2014 Return-path: Envelope-to: gllmg-musl@plane.gmane.org Original-Received: from mother.openwall.net ([195.42.179.200]) by plane.gmane.org with smtp (Exim 4.69) (envelope-from ) id 1XRyv0-0005XE-L5 for gllmg-musl@plane.gmane.org; Thu, 11 Sep 2014 09:35:50 +0200 Original-Received: (qmail 23962 invoked by uid 550); 11 Sep 2014 07:35:49 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: Original-Received: (qmail 23927 invoked from network); 11 Sep 2014 07:35:48 -0000 Content-Disposition: inline Xref: news.gmane.org gmane.linux.lib.musl.general:6135 Archived-At: --LZvS9be/3tNcYl/X Content-Type: text/plain; charset="us-ascii" Content-Disposition: inline Hello, Building musl with Clang against compiler-rt causes some of libc-test/math tests to fail on ARM. All failures are related to wrong floating point exceptions. Comparing GCC's builtins with their implementation in compiler-rt revealed important difference: compiler-rt builtins assume that input is always correct, while GCC builtins raise exceptions for invalid arguments or overflows/underflows. It looks like musl implicitly depends on GCC-like builtins. To fix that I added explicit checks in functions that failed to pass all tests. With these changes musl works in the same way independently whether it's built by Clang against compiler-rt or by GCC against libgcc. Please find the attached patch. It's only one as breaking it into pieces won't be very helpful. There is also separate patch for swapped arguments in error message of libc-test. Best regards, Sergey --LZvS9be/3tNcYl/X Content-Type: text/plain; charset="us-ascii" Content-Disposition: attachment; filename="0001-Raise-some-fenv-exceptions-explicitly.patch" >From cb8397f38a8a99884ffa85e2e4f6ec21e5d6cb87 Mon Sep 17 00:00:00 2001 From: Sergey Dmitrouk Date: Wed, 10 Sep 2014 15:09:49 +0300 Subject: [PATCH] Raise some fenv exceptions explicitly --- src/math/ilogb.c | 12 ++++++++++++ src/math/ilogbf.c | 12 ++++++++++++ src/math/j0.c | 9 +++++++-- src/math/j0f.c | 9 +++++++-- src/math/j1.c | 9 +++++++-- src/math/j1f.c | 9 +++++++-- src/math/jn.c | 5 ++++- src/math/jnf.c | 5 ++++- src/math/llrint.c | 13 ++++++++++++- src/math/llrintf.c | 13 ++++++++++++- src/math/llround.c | 13 ++++++++++++- src/math/llroundf.c | 13 ++++++++++++- src/math/llroundl.c | 13 ++++++++++++- src/math/pow.c | 33 +++++++++++++++++++++++++-------- src/math/sqrtl.c | 2 +- src/math/tgamma.c | 5 ++++- 16 files changed, 150 insertions(+), 25 deletions(-) diff --git a/src/math/ilogb.c b/src/math/ilogb.c index 64d4015..12c3aec 100644 --- a/src/math/ilogb.c +++ b/src/math/ilogb.c @@ -1,3 +1,5 @@ +#include +#include #include #include "libm.h" @@ -8,6 +10,16 @@ int ilogb(double x) uint64_t i = u.i; int e = i>>52 & 0x7ff; + if (x == 0.0) { + feraiseexcept(FE_INVALID); + return INT_MIN; + } + + if (isinf(x)) { + feraiseexcept(FE_INVALID); + return INT_MAX; + } + if (!e) { i <<= 12; if (i == 0) { diff --git a/src/math/ilogbf.c b/src/math/ilogbf.c index e23ba20..59fb2ba 100644 --- a/src/math/ilogbf.c +++ b/src/math/ilogbf.c @@ -1,4 +1,6 @@ +#include #include +#include #include "libm.h" int ilogbf(float x) @@ -8,6 +10,16 @@ int ilogbf(float x) uint32_t i = u.i; int e = i>>23 & 0xff; + if (x == 0.0) { + feraiseexcept(FE_INVALID); + return INT_MIN; + } + + if (isinf(x)) { + feraiseexcept(FE_INVALID); + return INT_MAX; + } + if (!e) { i <<= 9; if (i == 0) { diff --git a/src/math/j0.c b/src/math/j0.c index d722d94..9214d52 100644 --- a/src/math/j0.c +++ b/src/math/j0.c @@ -54,6 +54,7 @@ * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. */ +#include #include "libm.h" static double pzero(double), qzero(double); @@ -164,10 +165,14 @@ double y0(double x) EXTRACT_WORDS(ix, lx, x); /* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */ - if ((ix<<1 | lx) == 0) + if ((ix<<1 | lx) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0; + } if (ix >= 0x7ff00000) return 1/x; diff --git a/src/math/j0f.c b/src/math/j0f.c index 45883dc..2ac017b 100644 --- a/src/math/j0f.c +++ b/src/math/j0f.c @@ -14,6 +14,7 @@ */ #define _GNU_SOURCE +#include #include "libm.h" static float pzerof(float), qzerof(float); @@ -107,10 +108,14 @@ float y0f(float x) uint32_t ix; GET_FLOAT_WORD(ix, x); - if ((ix & 0x7fffffff) == 0) + if ((ix & 0x7fffffff) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0f; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0f; + } if (ix >= 0x7f800000) return 1/x; if (ix >= 0x40000000) { /* |x| >= 2.0 */ diff --git a/src/math/j1.c b/src/math/j1.c index df724d1..c83c568 100644 --- a/src/math/j1.c +++ b/src/math/j1.c @@ -54,6 +54,7 @@ * by method mentioned above. */ +#include #include "libm.h" static double pone(double), qone(double); @@ -156,10 +157,14 @@ double y1(double x) EXTRACT_WORDS(ix, lx, x); /* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */ - if ((ix<<1 | lx) == 0) + if ((ix<<1 | lx) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0; + } if (ix >= 0x7ff00000) return 1/x; diff --git a/src/math/j1f.c b/src/math/j1f.c index 58875af..3950d4d 100644 --- a/src/math/j1f.c +++ b/src/math/j1f.c @@ -14,6 +14,7 @@ */ #define _GNU_SOURCE +#include #include "libm.h" static float ponef(float), qonef(float); @@ -106,10 +107,14 @@ float y1f(float x) uint32_t ix; GET_FLOAT_WORD(ix, x); - if ((ix & 0x7fffffff) == 0) + if ((ix & 0x7fffffff) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0f; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0f; + } if (ix >= 0x7f800000) return 1/x; if (ix >= 0x40000000) /* |x| >= 2.0 */ diff --git a/src/math/jn.c b/src/math/jn.c index 4878a54..a9fa1cc 100644 --- a/src/math/jn.c +++ b/src/math/jn.c @@ -34,6 +34,7 @@ * values of n>1. */ +#include #include "libm.h" static const double invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */ @@ -224,8 +225,10 @@ double yn(int n, double x) if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */ return x; - if (sign && (ix|lx)!=0) /* x < 0 */ + if (sign && (ix|lx)!=0) { /* x < 0 */ + feraiseexcept(FE_INVALID); return 0/0.0; + } if (ix == 0x7ff00000) return 0.0; diff --git a/src/math/jnf.c b/src/math/jnf.c index f63c062..68b9cad 100644 --- a/src/math/jnf.c +++ b/src/math/jnf.c @@ -14,6 +14,7 @@ */ #define _GNU_SOURCE +#include #include "libm.h" float jnf(int n, float x) @@ -170,8 +171,10 @@ float ynf(int n, float x) ix &= 0x7fffffff; if (ix > 0x7f800000) /* nan */ return x; - if (sign && ix != 0) /* x < 0 */ + if (sign && ix != 0) { /* x < 0 */ + feraiseexcept(FE_INVALID); return 0/0.0f; + } if (ix == 0x7f800000) return 0.0f; diff --git a/src/math/llrint.c b/src/math/llrint.c index 4f583ae..701df94 100644 --- a/src/math/llrint.c +++ b/src/math/llrint.c @@ -1,8 +1,19 @@ +#include +#include #include /* uses LLONG_MAX > 2^53, see comments in lrint.c */ long long llrint(double x) { - return rint(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return x; + } + + x = rint(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llrintf.c b/src/math/llrintf.c index 96949a0..81153d2 100644 --- a/src/math/llrintf.c +++ b/src/math/llrintf.c @@ -1,8 +1,19 @@ +#include +#include #include /* uses LLONG_MAX > 2^24, see comments in lrint.c */ long long llrintf(float x) { - return rintf(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return x; + } + + x = rintf(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llround.c b/src/math/llround.c index 4d94787..ff48d9b 100644 --- a/src/math/llround.c +++ b/src/math/llround.c @@ -1,6 +1,17 @@ +#include +#include #include long long llround(double x) { - return round(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return 0; + } + + x = round(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llroundf.c b/src/math/llroundf.c index 19eb77e..b1bc9ec 100644 --- a/src/math/llroundf.c +++ b/src/math/llroundf.c @@ -1,6 +1,17 @@ +#include +#include #include long long llroundf(float x) { - return roundf(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return 0; + } + + x = roundf(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llroundl.c b/src/math/llroundl.c index 2c2ee5e..9a4e1a2 100644 --- a/src/math/llroundl.c +++ b/src/math/llroundl.c @@ -1,6 +1,17 @@ +#include +#include #include long long llroundl(long double x) { - return roundl(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return 0; + } + + x = roundl(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/pow.c b/src/math/pow.c index 82f684b..f19ea84 100644 --- a/src/math/pow.c +++ b/src/math/pow.c @@ -57,6 +57,7 @@ * to produce the hexadecimal values shown. */ +#include #include "libm.h" static const double @@ -196,16 +197,24 @@ double pow(double x, double y) /* |y| is huge */ if (iy > 0x41e00000) { /* if |y| > 2**31 */ if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */ - if (ix <= 0x3fefffff) + if (ix <= 0x3fefffff) { + feraiseexcept(FE_INEXACT| ((hy < 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy < 0 ? huge*huge : tiny*tiny; - if (ix >= 0x3ff00000) + } + if (ix >= 0x3ff00000) { + feraiseexcept(FE_INEXACT| ((hy > 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy > 0 ? huge*huge : tiny*tiny; + } } /* over/underflow if x is not close to one */ - if (ix < 0x3fefffff) + if (ix < 0x3fefffff) { + feraiseexcept(FE_INEXACT| ((hy < 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy < 0 ? s*huge*huge : s*tiny*tiny; - if (ix > 0x3ff00000) + } + if (ix > 0x3ff00000) { + feraiseexcept(FE_INEXACT| ((hy > 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy > 0 ? s*huge*huge : s*tiny*tiny; + } /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = ax - 1.0; /* t has 20 trailing zeros */ @@ -282,15 +291,23 @@ double pow(double x, double y) z = p_l + p_h; EXTRACT_WORDS(j, i, z); if (j >= 0x40900000) { /* z >= 1024 */ - if (((j-0x40900000)|i) != 0) /* if z > 1024 */ + if (((j-0x40900000)|i) != 0) { /* if z > 1024 */ + feraiseexcept(FE_INEXACT|FE_OVERFLOW); return s*huge*huge; /* overflow */ - if (p_l + ovt > z - p_h) + } + if (p_l + ovt > z - p_h) { + feraiseexcept(FE_INEXACT|FE_OVERFLOW); return s*huge*huge; /* overflow */ + } } else if ((j&0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */ // FIXME: instead of abs(j) use unsigned j - if (((j-0xc090cc00)|i) != 0) /* z < -1075 */ + if (((j-0xc090cc00)|i) != 0) { /* z < -1075 */ + feraiseexcept(FE_INEXACT|FE_UNDERFLOW); return s*tiny*tiny; /* underflow */ - if (p_l <= z - p_h) + } + if (p_l <= z - p_h) { + feraiseexcept(FE_INEXACT|FE_UNDERFLOW); return s*tiny*tiny; /* underflow */ + } } /* * compute 2**(p_h+p_l) diff --git a/src/math/sqrtl.c b/src/math/sqrtl.c index 83a8f80..0872e15 100644 --- a/src/math/sqrtl.c +++ b/src/math/sqrtl.c @@ -3,5 +3,5 @@ long double sqrtl(long double x) { /* FIXME: implement in C, this is for LDBL_MANT_DIG == 64 only */ - return sqrt(x); + return isnan(x) ? 0.0l/0.0l : sqrt(x); } diff --git a/src/math/tgamma.c b/src/math/tgamma.c index 28f6e0f..51e7e20 100644 --- a/src/math/tgamma.c +++ b/src/math/tgamma.c @@ -22,6 +22,7 @@ Gamma(x)*Gamma(-x) = -pi/(x sin(pi x)) most ideas and constants are from boost and python */ +#include #include "libm.h" static const double pi = 3.141592653589793238462643383279502884; @@ -124,8 +125,10 @@ double tgamma(double x) /* integer arguments */ /* raise inexact when non-integer */ if (x == floor(x)) { - if (sign) + if (sign) { + feraiseexcept(FE_INVALID); return 0/0.0; + } if (x <= sizeof fact/sizeof *fact) return fact[(int)x - 1]; } -- 1.8.4 --LZvS9be/3tNcYl/X Content-Type: text/plain; charset="us-ascii" Content-Disposition: attachment; filename="0001-Fix-order-of-jn-arguments-in-error-message.patch" >From 7b5026b43144faff415c948f47df40c3a9d5dc81 Mon Sep 17 00:00:00 2001 From: Sergey Dmitrouk Date: Wed, 10 Sep 2014 14:33:42 +0300 Subject: [PATCH] Fix order of jn() arguments in error message They are swapped. --- src/math/jn.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/math/jn.c b/src/math/jn.c index 7573edc..e02db5e 100644 --- a/src/math/jn.c +++ b/src/math/jn.c @@ -27,17 +27,17 @@ int main(void) y = jn(p->i, p->x); e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW); if (!checkexcept(e, p->e, p->r)) { - printf("%s:%d: bad fp exception: %s jn(%a, %lld)=%a, want %s", - p->file, p->line, rstr(p->r), p->x, p->i, p->y, estr(p->e)); + printf("%s:%d: bad fp exception: %s jn(%lld, %a)=%a, want %s", + p->file, p->line, rstr(p->r), p->i, p->x, p->y, estr(p->e)); printf(" got %s\n", estr(e)); err++; } d = ulperr(y, p->y, p->dy); if (!checkulp(d, p->r)) { - printf("%s:%d: %s jn(%a, %lld) want %a got %a, ulperr %.3f = %a + %a\n", - p->file, p->line, rstr(p->r), p->x, p->i, p->y, y, d, d-p->dy, p->dy); + printf("%s:%d: %s jn(%lld, %a) want %a got %a, ulperr %.3f = %a + %a\n", + p->file, p->line, rstr(p->r), p->i, p->x, p->y, y, d, d-p->dy, p->dy); err++; } } return !!err; -- 1.8.4 --LZvS9be/3tNcYl/X--