[PATCH] Make musl math depend less on libgcc builtins

[Sergey Dmitrouk <sdmitrouk@accesssoftek.com>]

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