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

[Szabolcs Nagy <nsz@port70.net>]

* 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