[PATCH] Make musl math depend less on libgcc builtins

[PATCH] Make musl math depend less on libgcc builtins

[Sergey Dmitrouk]

[-- Attachment #1: Type: text/plain, Size: 886 bytes --]

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

[-- Attachment #2: 0001-Raise-some-fenv-exceptions-explicitly.patch --]
[-- Type: text/plain, Size: 11105 bytes --]

From cb8397f38a8a99884ffa85e2e4f6ec21e5d6cb87 Mon Sep 17 00:00:00 2001
From: Sergey Dmitrouk <sdmitrouk@accesssoftek.com>
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 <fenv.h>
+#include <math.h>
 #include <limits.h>
 #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 <fenv.h>
 #include <limits.h>
+#include <math.h>
 #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 <fenv.h>
 #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 <fenv.h>
 #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 <fenv.h>
 #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 <fenv.h>
 #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 <fenv.h>
 #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 <fenv.h>
 #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 <fenv.h>
+#include <limits.h>
 #include <math.h>
 
 /* 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 <fenv.h>
+#include <limits.h>
 #include <math.h>
 
 /* 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 <fenv.h>
+#include <limits.h>
 #include <math.h>
 
 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 <fenv.h>
+#include <limits.h>
 #include <math.h>
 
 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 <fenv.h>
+#include <limits.h>
 #include <math.h>
 
 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 <fenv.h>
 #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 <fenv.h>
 #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


[-- Attachment #3: 0001-Fix-order-of-jn-arguments-in-error-message.patch --]
[-- Type: text/plain, Size: 1330 bytes --]

From 7b5026b43144faff415c948f47df40c3a9d5dc81 Mon Sep 17 00:00:00 2001
From: Sergey Dmitrouk <sdmitrouk@accesssoftek.com>
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

Re: [PATCH] Make musl math depend less on libgcc builtins

[Szabolcs Nagy]

* Sergey Dmitrouk <sdmitrouk@accesssoftek.com> [2014-09-11 10:35:32 +0300]:
> 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.

floating-point arithmetics is not done by the compiler runtime

(unless it is software emulated, but then exceptions are not implemented
by libgcc because it used to be impossible without hw registers for it
in c99, the signal handler semantics changed in c11 and now fenv is
thread local storage so they could implement exceptions through tls...)

the exception semantics is part of the ieee 754 spec (which is referenced
by iso c annex f) so there cannot be a difference between gcc and clang
(if there is then one of them is broken.. or both of them)

it turns out neither gcc nor clang imlpement annex f correctly (there are
famous bug reports for it on both sides), in particular they misoptimize
code that depends on non-default fenv or accesses fenv

gcc just happens to work usually with -frounding-math, clang makes no attempt
to undo its optimizations with any compiler flags (except -O0 which is not
what most ppl want) like constant folding 1/0.0 into INFINITY

this is a sad issue (just ask freebsd math maintainers, i think they had
to hunt down a lot of miscompilations in lib/msun to make clang work.
in musl we rather not support broken compilers by littering the code
with volatile hacks or asm memory barriers to get reasonable behaviour)

> 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.

i dont see any gcc builtins used, we rely on c99 fenv semantics

if a compiler does not support

 #pragma STDC FENV_ACCESS ON

then it has to be conservative about optimizations and assume that
fenv access is always on.

(neither gcc nor clang supports FENV_ACCESS, gcc just happens to
be close to always on semantics with -frounding-math)

> 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

> >From cb8397f38a8a99884ffa85e2e4f6ec21e5d6cb87 Mon Sep 17 00:00:00 2001
> From: Sergey Dmitrouk <sdmitrouk@accesssoftek.com>
> 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 <fenv.h>
> +#include <math.h>
>  #include <limits.h>
>  #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;
> +	}
> +

well this has some issues

FE_* macros may be undefined for a target so their use always
have to be ifdefed

using feraiseexcept is slower and bigger than relying on the
fpu doing the exception raising for us (this may not be a big
problem since we only support inexact exception where absolutely
necessary and other exceptions usually occur in rare code paths)

feraiseexcept also adds extra dependency so math functions
cannot be used without a working fenv then

the compiler can still miscompile the code: feraiseexcept may
be implemented using arithmetics and then we have the same
issue

sometimes it's not clear if the compiler can optimize some
arithmetics, eg:

 y = 0*x;

may be turned into

 y = isinf(x) ? NAN : 0;

so should we raise the invalid flag manually or rely on that
the compiler will use fpu intructions which do it for us?

that said, i'm open to changes in the current policy since no
compiler supports FENV_ACCESS correctly and there does not seem
to be much willingness to fix this

- reorderings around fesetround, fetestexcept, feclearexcept
are harder to fix, but we only use those in a few places so
volatile hacks may not be terrible

- for exception raising if we can reliably identify the places
where the compiler miscompiles/constant folds the code then we
can fix those with explicit feraise (or volatile hacks) if it
does not have too much impact otherwise

> 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 <fenv.h>
> +#include <limits.h>
>  #include <math.h>
>  
>  /* 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;
>  }

this should not be needed, overflowing float to int conversion
raises the invalid flag implicitly, if it does not then clang/llvm
generates wrong code for the conversion

> 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);
>  }

why?

0/.0 raises invalid exception but quiet nan never raises
exception (unless converted to int or used in <,> relations)

nan is also sticky (passes through any arithmetics and
comes out as nan) so if sqrt(NAN) is not nan now then
that's a bug somewhere

> >From 7b5026b43144faff415c948f47df40c3a9d5dc81 Mon Sep 17 00:00:00 2001
> From: Sergey Dmitrouk <sdmitrouk@accesssoftek.com>
> Date: Wed, 10 Sep 2014 14:33:42 +0300
> Subject: [PATCH] Fix order of jn() arguments in error message
> 
> They are swapped.
> ---

applied and did the same for jnf, yn, ynf

Re: [PATCH] Make musl math depend less on libgcc builtins

[Sergey Dmitrouk]

Hi,

> floating-point arithmetics is not done by the compiler runtime

I do understand that, but got lost after observing different results
produced by same code compiled with different compilers.

> misoptimize code that depends on non-default fenv or accesses fenv
>
> gcc just happens to work usually with -frounding-math, clang makes no attempt
> to undo its optimizations with any compiler flags (except -O0 which is not
> what most ppl want) like constant folding 1/0.0 into INFINITY

I though that lack of -frounding-math flag support in Clang might be the
reason of errors, but I didn't realize optimizations can harm floating
point operations.

> i dont see any gcc builtins used, we rely on c99 fenv semantics

I meant builtins inserted by compiler, which I saw while debugging
tests.  Now I think that was soft floating point, although it's not what
I want, hard float should be used instead.  I might have relied to much
on compiler defaults.

> FE_* macros may be undefined for a target so their use always
> have to be ifdefed

I read this somewhere, but wasn't sure that it's target specific.

> so should we raise the invalid flag manually or rely on that
> the compiler will use fpu intructions which do it for us?
> 
> that said, i'm open to changes in the current policy since no
> compiler supports FENV_ACCESS correctly and there does not seem
> to be much willingness to fix this
> 
> - reorderings around fesetround, fetestexcept, feclearexcept
> are harder to fix, but we only use those in a few places so
> volatile hacks may not be terrible
> 
> - for exception raising if we can reliably identify the places
> where the compiler miscompiles/constant folds the code then we
> can fix those with explicit feraise (or volatile hacks) if it
> does not have too much impact otherwise

Fixing it in musl won't help other applications compiled with Clang, so
I'd prefer to fix such issues in the compiler.

> this should not be needed, overflowing float to int conversion
> raises the invalid flag implicitly, if it does not then clang/llvm
> generates wrong code for the conversion

Well, it doesn't, will need to figure out why.  Because of strange
results I interpreted the whole thing in a wrong way, as if exceptions
were semi automatic and libc implementation had to raise some exceptions
manually in places where things defined by C standard differ from what
IEEE754 implementation gives us.  You helpful comments sorted that out
for me.

> > 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);
> >  }
>
> why?

sqrt(NAN) raises INVALID exception, 0.0l/0.0l doesn't for me (well,
optimization must've prevented that).

> nan is also sticky (passes through any arithmetics and
> comes out as nan) so if sqrt(NAN) is not nan now then
> that's a bug somewhere

sqrt(NAN) == NAN, I just wanted to silent the exception.

> applied and did the same for jnf, yn, ynf

Thanks, I didn't notice it in other tests.


Thanks for your verbose comments.  I'll need to figure out why exceptions
are not raised in some cases, but now I understand that I went into wrong
direction blaming musl on this.

Regards,
Sergey

Re: [PATCH] Make musl math depend less on libgcc builtins

[Szabolcs Nagy]

* Sergey Dmitrouk <sdmitrouk@accesssoftek.com> [2014-09-11 14:42:07 +0300]:
> > gcc just happens to work usually with -frounding-math, clang makes no attempt
> > to undo its optimizations with any compiler flags (except -O0 which is not
> > what most ppl want) like constant folding 1/0.0 into INFINITY
> 
> I though that lack of -frounding-math flag support in Clang might be the
> reason of errors, but I didn't realize optimizations can harm floating
> point operations.
> 

note that constant folding 1/0.0 is valid if FENV_ACCESS is OFF
(it is off by default), however if FENV_ACCESS ON is not supported
then ON should be always assumed and then constant folding is not ok
(this is where both gcc and clang are wrong when -std=c99 is specified
however they don't claim complete annex f support so it is hard to
argue this point)

http://port70.net/~nsz/c/c11/n1570.html#F.8.4p1

> > this should not be needed, overflowing float to int conversion
> > raises the invalid flag implicitly, if it does not then clang/llvm
> > generates wrong code for the conversion
> 
> Well, it doesn't, will need to figure out why.  Because of strange
> results I interpreted the whole thing in a wrong way, as if exceptions
> were semi automatic and libc implementation had to raise some exceptions
> manually in places where things defined by C standard differ from what
> IEEE754 implementation gives us.  You helpful comments sorted that out
> for me.
> 
> > > 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);
> > >  }
> >
> > why?
> 
> sqrt(NAN) raises INVALID exception, 0.0l/0.0l doesn't for me (well,
> optimization must've prevented that).
> 
> > nan is also sticky (passes through any arithmetics and
> > comes out as nan) so if sqrt(NAN) is not nan now then
> > that's a bug somewhere
> 
> sqrt(NAN) == NAN, I just wanted to silent the exception.

now that you mention i got reports of spurious invalid
exceptions for nan inputs on arm, but i don't have hardfloat
arm toolchain

if you have further info on this that would be helpful

(from the arm docs it seems to me that vfp sqrt should work
according to ieee so either that is wrong or our fenv
exception test is not ok on arm)

Re: [PATCH] Make musl math depend less on libgcc builtins

[Sergey Dmitrouk]

> now that you mention i got reports of spurious invalid
> exceptions for nan inputs on arm, but i don't have hardfloat
> arm toolchain
> 
> if you have further info on this that would be helpful
> 
> (from the arm docs it seems to me that vfp sqrt should work
> according to ieee so either that is wrong or our fenv
> exception test is not ok on arm)

Note that I'm performing tests in QEMU.  vsqrt.f64 instruction returns
NAN for NAN input and doesn't raise any exceptions ($fpscr register is
unchanged).  However I get INVALID exception after printing result with
musl implementation of printf(), glibc doesn't raise any exceptions in
printf(), but I'm not sure whether it's important, probably not.

Re: [PATCH] Make musl math depend less on libgcc builtins

[Szabolcs Nagy]

* Sergey Dmitrouk <sdmitrouk@accesssoftek.com> [2014-09-11 16:22:53 +0300]:
> Note that I'm performing tests in QEMU.  vsqrt.f64 instruction returns
> NAN for NAN input and doesn't raise any exceptions ($fpscr register is
> unchanged).  However I get INVALID exception after printing result with
> musl implementation of printf(), glibc doesn't raise any exceptions in
> printf(), but I'm not sure whether it's important, probably not.

nan printf should not raise invalid exceptino in musl either

musl checks for nan by if(y!=y) and my guess is that this is
incorrectly done by a signaling comparision

i checked on x86_64 and both clang and gcc get comparisions
wrong in the other direction: they use quite comparisions
when signaling is needed, eg

http://goo.gl/GsdpZA

(on a correct implementation ==, != are quiet, but <,>,<=,>=
raise invalid if any of the operands are nan, on x86_64 the
quiet instruction is ucomis* and the signaling one is comis*
and both gcc-4.9 and clang-3.4 seem to use ucomis* for all
relational operations, may be there is some compiler flag to
make them behave..)

Re: [PATCH] Make musl math depend less on libgcc builtins

[Sergey Dmitrouk]

On Thu, Sep 11, 2014 at 07:11:23AM -0700, Szabolcs Nagy wrote:
> musl checks for nan by if(y!=y) and my guess is that this is
> incorrectly done by a signaling comparision

I see

    vcmpe.f64      d0, d0

which seems to be part of "!isfinite(y)" expression, so your guess is
correct ("vcmpe" instruction raises exceptions, "vcmp" doesn't).

I also checked >, <, <= and >=, which Clang implements as "vcmpe" too.
Linaro GCC 4.9 seems to do the correct thing for ARMhf.  So different
targets behave differently, maybe there is no flag to control this...

> 
> i checked on x86_64 and both clang and gcc get comparisions
> wrong in the other direction: they use quite comparisions
> when signaling is needed, eg
> 
> http://goo.gl/GsdpZA
> 
> (on a correct implementation ==, != are quiet, but <,>,<=,>=
> raise invalid if any of the operands are nan, on x86_64 the
> quiet instruction is ucomis* and the signaling one is comis*
> and both gcc-4.9 and clang-3.4 seem to use ucomis* for all
> relational operations, may be there is some compiler flag to
> make them behave..)

Re: [PATCH] Make musl math depend less on libgcc builtins

[Alexander Monakov]

On Thu, 11 Sep 2014, Szabolcs Nagy wrote:

> * Sergey Dmitrouk <sdmitrouk@accesssoftek.com> [2014-09-11 16:22:53 +0300]:
> > Note that I'm performing tests in QEMU.  vsqrt.f64 instruction returns
> > NAN for NAN input and doesn't raise any exceptions ($fpscr register is
> > unchanged).  However I get INVALID exception after printing result with
> > musl implementation of printf(), glibc doesn't raise any exceptions in
> > printf(), but I'm not sure whether it's important, probably not.
> 
> nan printf should not raise invalid exceptino in musl either
> 
> musl checks for nan by if(y!=y) and my guess is that this is
> incorrectly done by a signaling comparision
> 
> i checked on x86_64 and both clang and gcc get comparisions
> wrong in the other direction: they use quite comparisions
> when signaling is needed, eg
> 
> http://goo.gl/GsdpZA
> 
> (on a correct implementation ==, != are quiet, but <,>,<=,>=
> raise invalid if any of the operands are nan, on x86_64 the
> quiet instruction is ucomis* and the signaling one is comis*
> and both gcc-4.9 and clang-3.4 seem to use ucomis* for all
> relational operations, may be there is some compiler flag to
> make them behave..)

In GCC, -mno-ieee-fp will switch ucom->com it for both == and <>.
Is that still wrong?  Shouldn't there be a GCC bug filed?

Alexander

Re: [PATCH] Make musl math depend less on libgcc builtins

[Szabolcs Nagy]

* Alexander Monakov <amonakov@ispras.ru> [2014-09-11 19:14:50 +0400]:
> 
> In GCC, -mno-ieee-fp will switch ucom->com it for both == and <>.
> Is that still wrong?  Shouldn't there be a GCC bug filed?
> 

yes

i's backwards for <> ("ieee-fp" should mean ieee semantics
and "no-ieee-fp" something else) and ==, != needs to be ucom

Re: [PATCH] Make musl math depend less on libgcc builtins

[Sergey Dmitrouk]

On Thu, Sep 11, 2014 at 07:11:23AM -0700, Szabolcs Nagy wrote:
> on a correct implementation ==, != are quiet, but <,>,<=,>=
> raise invalid if any of the operands are nan

I wanted to get some details on this, but failed to find relevant
sections of C99/IEEE754 standards.  I see C99 referring to IEEE754, is
it in "5.11 Details of comparison predicates" section of IEEE754?
Could you please point me to section(s) I'm apparently missing?

-- 
Sergey

Re: [PATCH] Make musl math depend less on libgcc builtins

[Szabolcs Nagy]

* Sergey Dmitrouk <sdmitrouk@accesssoftek.com> [2014-09-18 17:28:40 +0300]:
> On Thu, Sep 11, 2014 at 07:11:23AM -0700, Szabolcs Nagy wrote:
> > on a correct implementation ==, != are quiet, but <,>,<=,>=
> > raise invalid if any of the operands are nan
> 
> I wanted to get some details on this, but failed to find relevant
> sections of C99/IEEE754 standards.  I see C99 referring to IEEE754, is
> it in "5.11 Details of comparison predicates" section of IEEE754?
> Could you please point me to section(s) I'm apparently missing?

yes ieee754-2008 5.11

i think it is clear: there are tables showing all the predicates
and to which "traditional names and symbols" they should map.

table 5.1 shows ==, != as quiet comparisions, table 5.2 shows
<,> operations as signaling and the text mentions that the quiet
operations in table 5.3 are for applications which want to
explicitly handle quiet nans that way

the text in iso C F.3 is not very detailed about the mapping but
gives hints:

 The relational and equality operators provide IEC 60559 comparisons.
 IEC 60559 identifies a need for additional comparison predicates to
 facilitate writing code that accounts for NaNs. The comparison macros
 (isgreater, isgreaterequal, isless, islessequal, islessgreater, and
 isunordered) in <math.h> supplement the language operators to address
 this need. The islessgreater and isunordered macros provide respectively
 a quiet version of the <> predicate and the unordered predicate
 recommended in the Appendix to IEC 60559.

the <,> predicates need a quiet version because the default is not quiet,
but == and != dont since they are already quiet

the precise mapping will be spelled out in more detail in TS 18661,
 http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1778.pdf
see "table - 1 operation binding"
(the latest version seems to be password protected, sigh..)

Re: [PATCH] Make musl math depend less on libgcc builtins

[Sergey Dmitrouk]

On Thu, Sep 18, 2014 at 10:55:47AM -0700, Szabolcs Nagy wrote:
> i think it is clear: there are tables showing all the predicates
> and to which "traditional names and symbols" they should map.

"EQ" predicate appearing in both 5.1 and 5.2 tables confused me, I
missed that "=" in 5.1 and no "=" in 5.2 means binding to actual
operation.

> table 5.1 shows ==, != as quiet comparisions, table 5.2 shows
> <,> operations as signaling and the text mentions that the quiet
> operations in table 5.3 are for applications which want to
> explicitly handle quiet nans that way
> 
> the text in iso C F.3 is not very detailed about the mapping but
> gives hints:
> 
>  The relational and equality operators provide IEC 60559 comparisons.
>  IEC 60559 identifies a need for additional comparison predicates to
>  facilitate writing code that accounts for NaNs. The comparison macros
>  (isgreater, isgreaterequal, isless, islessequal, islessgreater, and
>  isunordered) in <math.h> supplement the language operators to address
>  this need. The islessgreater and isunordered macros provide respectively
>  a quiet version of the <> predicate and the unordered predicate
>  recommended in the Appendix to IEC 60559.
> 
> the <,> predicates need a quiet version because the default is not quiet,
> but == and != dont since they are already quiet
> 
> the precise mapping will be spelled out in more detail in TS 18661,
>  http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1778.pdf
> see "table - 1 operation binding"
> (the latest version seems to be password protected, sigh..)

Thanks for the explanation.