double f(){ return __builtin_exp(-__builtin_huge_val()); } As noticed in PR57974, we fail to fold this to a constant. Indeed, in do_mpfr_arg1, the relevant code is protected by if (real_isfinite (ra) ... mpfr handles infinite values, so it should be doable, at least with -ffast-math (there might be some errno issues with default options).

Thanks Marc.

There are no errno issues - this is an exact zero result, not underflow. But I'm not confident that MPFR follows all the Annex F special cases for infinities and NaNs (and even less confident in MPC following Annex G), and in cases that are errors and should produce exceptions / errno (e.g. sin (Inf)) you do of course need to avoid folding.

Oh nice. And if I disable by hand the real_isfinite (ra) check in do_mpfr_arg1 I even get 0. And I also checked what happens for sin(Inf) in that case: a -nan as before the hack. Then which is at this point a safe way to proceed? Get in touch with the mpfr people, ask if simplifying infinities has known issues? Tentatively remove the real_isfinite checks from one of the do_mpfr_arg? functions at a time, or even one mathematical function at a time?

The MPFR documentation does claim that it strictly conforms to annex F (with an explanation on how to emulate subnormals), though it isn't clear if that claim only concerns +-*/sqrt or everything. (In reply to joseph@codesourcery.com from comment #2) > There are no errno issues - this is an exact zero result, not underflow. > But I'm not confident that MPFR follows all the Annex F special cases for > infinities and NaNs (and even less confident in MPC following Annex G), > and in cases that are errors and should produce exceptions / errno (e.g. > sin (Inf)) you do of course need to avoid folding. I was wondering about that last point. Couldn't we replace: x=sin(Inf); with: x=NaN; errno=EDOM; // only if flag_math_errno volatile double f=NaN+NaN; // if flag_trapping_math, something to raise invalid (make sure we don't recursively try to propagate the constant there, so maybe the NaN argument should be volatile) Ok, NaN may not be the most interesting case to propagate, but infinities are interesting, even when they would give EDOM / FE_OVERFLOW. Though the compiler might need help from the front-end finding the variable errno and the values EDOM and FE_OVERFLOW if it needs them explicitly. Of course, doing the infinity propagation only in the !flag_math_errno && !flag_trapping_math (and I guess !flag_rounding_math) case is a natural first step.

Today I was thinking that given that, per docs and testsuite (double checked yesterday) the mpfr functions are able to cope with +-Inf arguments to the mathematical functions and evaluate correctly, gating the various do_mpfr_* with !real_isnan instead of real_isfinite doesn't look like taking a big risk, now that we are in Stage 1. Alone that would help a lot of code (in particular, in the C++ library, which is the original motivating example). Note, I'm not thinking replacing real_isfinite in various other places, in particular not in do_mpc_* and its helpers (something for the future). Comments?

(In reply to Paolo Carlini from comment #5) > Today I was thinking that given that, per docs and testsuite (double checked > yesterday) the mpfr functions are able to cope with +-Inf arguments to the > mathematical functions and evaluate correctly, gating the various do_mpfr_* > with !real_isnan instead of real_isfinite doesn't look like taking a big > risk, now that we are in Stage 1. Alone that would help a lot of code (in > particular, in the C++ library, which is the original motivating example). > Note, I'm not thinking replacing real_isfinite in various other places, in > particular not in do_mpc_* and its helpers (something for the future). > Comments? I haven't looked at that code closely enough to say if there is much that needs changing. Earlier you said that removing the real_isfinite test replaced sin(Inf) with NaN. Was that with -ffast-math? Without it, the change shouldn't happen (not without inserting a few extra statements). I am hoping that there is already code in place that checks if the mpfr calls set errno and calls mpfr_overflow_p (and family) for possible exceptions, though I somehow doubt that we test mpfr_inexflag_p or very little propagation would take place.

An example of MPC not following all the Annex G special cases is that catanh (1 + i0) is specified in Annex G to return Inf + i0 with divide-by-zero exception, but at least with my MPC installation it returns Inf + iNaN. I haven't tried to check how MPFR handles special cases; the issue with MPC is simply something I noticed incidentally while fixing glibc libm handling of various <complex.h> functions. > I was wondering about that last point. Couldn't we replace: > > x=sin(Inf); > > with: > > x=NaN; > errno=EDOM; // only if flag_math_errno errno is typically a macro.... > volatile double f=NaN+NaN; // if flag_trapping_math, something to raise invalid > (make sure we don't recursively try to propagate the constant there, so maybe > the NaN argument should be volatile) I think you mean 0.0 / 0.0 or Inf - Inf or similar; NaN+NaN doesn't raise an exception if the NaNs are quiet NaNs.

About sin(Inf): I checked that with / without the real_isfinite the expression evaluated in every case the same, -nan, if I remember correctly. I don't have more details at the moment. My general point, again, is that the mpfr testsuite appears to include test which check that +-inf are correctly handled as arguments to the mathematical functions. Which, hey, doesn't seem a miracle, after all ;) Remember, I'm thinking as a start, no-nan, no-complex.

(In reply to joseph@codesourcery.com from comment #7) > An example of MPC not following all the Annex G special cases is that > catanh (1 + i0) is specified in Annex G to return Inf + i0 with > divide-by-zero exception, but at least with my MPC installation it returns > Inf + iNaN. I haven't tried to check how MPFR handles special cases; the > issue with MPC is simply something I noticed incidentally while fixing > glibc libm handling of various <complex.h> functions. Thanks (I assume you reported it to MPC, so that will be one fewer issue in a few years :-). I believe there are far fewer special cases (and thus risks) with MPFR, but that would indeed require a suitable testsuite for all functions for which we enable it (at least if MPFR doesn't already have such a testsuite, and maybe even then, to make sure we call it properly). > > I was wondering about that last point. Couldn't we replace: > > > > x=sin(Inf); > > > > with: > > > > x=NaN; > > errno=EDOM; // only if flag_math_errno > > errno is typically a macro.... That's why I was mentioning front-end help... There should be ways to set errno to EDOM faster than calling sin(Inf). > > volatile double f=NaN+NaN; // if flag_trapping_math, something to raise invalid > > (make sure we don't recursively try to propagate the constant there, so maybe > > the NaN argument should be volatile) > > I think you mean 0.0 / 0.0 or Inf - Inf or similar; NaN+NaN doesn't raise > an exception if the NaNs are quiet NaNs. Yeah, any of those. I was inspired by glibc, which has for instance: double __fdim (double x, double y) { int clsx = fpclassify (x); int clsy = fpclassify (y); if (clsx == FP_NAN || clsy == FP_NAN || (y < 0 && clsx == FP_INFINITE && clsy == FP_INFINITE)) /* Raise invalid flag. */ return x - y; which looks like it expects QNaN-QNaN to set the invalid flag.

About testing, it would be just matter of extending/updating what Kaveh Ghazi set up when mpfr/mpc came in.

On Sat, 27 Jul 2013, glisse at gcc dot gnu.org wrote: > Yeah, any of those. I was inspired by glibc, which has for instance: > > double > __fdim (double x, double y) > { > int clsx = fpclassify (x); > int clsy = fpclassify (y); > > if (clsx == FP_NAN || clsy == FP_NAN > || (y < 0 && clsx == FP_INFINITE && clsy == FP_INFINITE)) > /* Raise invalid flag. */ > return x - y; > > which looks like it expects QNaN-QNaN to set the invalid flag. Such comments must be understood to be written on the assumption that the reader is familiar with the desired IEEE semantics - that is, that the flags is raised if and only if a NaN argument is a signaling NaN (and such arithmetic patterns on input NaNs, to ensure that "invalid" is raised if either NaN is signaling, and otherwise that an input NaN's significand is preserved, are very common in glibc).

(In reply to Marc Glisse from comment #9) > I believe there are far fewer special cases (and thus > risks) with MPFR, but that would indeed require a suitable testsuite for all > functions for which we enable it (at least if MPFR doesn't already have such > a testsuite, and maybe even then, to make sure we call it properly). MPFR's testsuite is just against the MPFR implementation. These are actually non-regression tests. For comparisons with the functions from the C library, there's mpcheck: https://gforge.inria.fr/projects/mpcheck/ but I don't know whether it includes special values (it wasn't its goal). Note that we tried to follow C99's Annex F when this made sense. MPFR also supports some special functions that are not in ISO C (yet), but may be provided by the C library on some platforms (e.g. Bessel functions, which are also specified in POSIX). Don't forget that the specific rules for signed zeroes are also concerned; again, we tried to follow C99's Annex F, IIRC, even when the specification was rather inconsistent (e.g. under the rounding mode toward negative infinity, the subtraction 1 - 1 returns -0, but log(1) is required to return +0).

A difference that may occur in the future if the C library adds a rsqrt function (based on the IEEE 754-2008 rSqrt function) or constant folding is used on builtins: in MPFR, mpfr_rec_sqrt on -0 gives +Inf, not -Inf.

Let's see.

If GCC intends to handle Bessel functions j0, j1, jn, y0, y1, yn (POSIX), there may be differences with GNU MPFR. See my messages and bug report: http://permalink.gmane.org/gmane.comp.standards.posix.austin.general/7982 http://permalink.gmane.org/gmane.comp.standards.posix.austin.general/7990 http://sourceware.org/bugzilla/show_bug.cgi?id=15901

Uhm, something isn't going as planned. If I change that check in do_mpfr_arg1 to !real_isnan (ra), then we always fold __builtin_expl(-__builtin_huge_vall()) to an exact 0. However, with *both* -fno-math-errno -funsafe-math-optimizations we don't fold the kind of code we really care about: type num = __builtin_logl(0.0L); type res = __builtin_expl(num); and res turns out to be -nan :( Probably the problem has to do with the log, not sure yet.

If I play with some constexprs, I understand that, irrespective of the command line switches, we never fold to a constant the log. However, only for -fno-math-errno -funsafe-math-optimizations the library call returns -nan, instead of 0. Ok.. For the log, the check at the beginning of do_mpfr_arg1: if (!real_isnan (ra) && (!min || real_compare (inclusive ? GE_EXPR: GT_EXPR , ra, min)) && (!max || real_compare (inclusive ? LE_EXPR: LT_EXPR , ra, max))) is false, because real_compare (inclusive ? GE_EXPR: GT_EXPR , ra, min) is false. It seems we never fold log(0).

Now the question is: why we use a library call for log(0) instead of mpfr?

If I change fold_builtin_logarithm to pass a true as last argument to do_mpfr_arg1 (thus 0 is accepted) and do_mpfr_ckconv to accept a folded result which is infinity, things finally work. Patchlet below. Note however, that I also need -O1 otherwise, at -O0, we don't try to propagate the constant num and mpfr isn't used, we again have a library call which returns -nan. Index: builtins.c =================================================================== --- builtins.c (revision 204005) +++ builtins.c (working copy) @@ -8191,7 +8191,7 @@ fold_builtin_logarithm (location_t loc, tree fndec const enum built_in_function fcode = builtin_mathfn_code (arg); /* Calculate the result when the argument is a constant. */ - if ((res = do_mpfr_arg1 (arg, type, func, &dconst0, NULL, false))) + if ((res = do_mpfr_arg1 (arg, type, func, &dconst0, NULL, true))) return res; /* Special case, optimize logN(expN(x)) = x. */ @@ -13527,7 +13527,7 @@ do_mpfr_ckconv (mpfr_srcptr m, tree type, int inex /* Proceed iff we get a normal number, i.e. not NaN or Inf and no overflow/underflow occurred. If -frounding-math, proceed iff the result of calling FUNC was exact. */ - if (mpfr_number_p (m) && !mpfr_overflow_p () && !mpfr_underflow_p () + if (!mpfr_nan_p (m) && !mpfr_overflow_p () && !mpfr_underflow_p () && (!flag_rounding_math || !inexact)) { REAL_VALUE_TYPE rr; @@ -13537,7 +13537,7 @@ do_mpfr_ckconv (mpfr_srcptr m, tree type, int inex check for overflow/underflow. If the REAL_VALUE_TYPE is zero but the mpft_t is not, then we underflowed in the conversion. */ - if (real_isfinite (&rr) + if (!real_isnan (&rr) && (rr.cl == rvc_zero) == (mpfr_zero_p (m) != 0)) { REAL_VALUE_TYPE rmode; @@ -13623,7 +13623,7 @@ do_mpfr_arg1 (tree arg, tree type, int (*func)(mpf { const REAL_VALUE_TYPE *const ra = &TREE_REAL_CST (arg); - if (real_isfinite (ra) + if (!real_isnan (ra) && (!min || real_compare (inclusive ? GE_EXPR: GT_EXPR , ra, min)) && (!max || real_compare (inclusive ? LE_EXPR: LT_EXPR , ra, max))) {

Thus I clearly realize something: if we do better for mpfr, the issue remains that if we do *not* optimize and constants are not propagated, we issue library calls, which, for command line switches like -fno-math-errno -funsafe-math-optimizations give incorrect nan result. In other terms, I'm still more than willing to work on mpfr-related improvements, but I'm still not convinced that when constants are not at issue we are really ok with those nans... Note again, in case isn't clear already, that for this kind of snippet I never managed to see nans with clang/icc. About mpfr, I would appreciate feedback that I'm on the right track. Thanks!

(In reply to Paolo Carlini from comment #19) > If I change fold_builtin_logarithm to pass a true as last argument to > do_mpfr_arg1 (thus 0 is accepted) and do_mpfr_ckconv to accept a folded > result which is infinity, things finally work. Patchlet below. Note however, > that I also need -O1 otherwise, at -O0, we don't try to propagate the > constant num and mpfr isn't used, we again have a library call which returns > -nan. I don't know much about GCC internals, but did you think about exceptions, e.g. FE_DIVBYZERO for log(0)? Note that MPFR's divide-by-zero exception (flag) and associated functions are new in MPFR 3.1.0. > @@ -13527,7 +13527,7 @@ do_mpfr_ckconv (mpfr_srcptr m, tree type, int inex > /* Proceed iff we get a normal number, i.e. not NaN or Inf and no > overflow/underflow occurred. If -frounding-math, proceed iff the > result of calling FUNC was exact. */ > - if (mpfr_number_p (m) && !mpfr_overflow_p () && !mpfr_underflow_p () > + if (!mpfr_nan_p (m) && !mpfr_overflow_p () && !mpfr_underflow_p () > && (!flag_rounding_math || !inexact)) [...] If you do this, don't forget to update the comment. Ditto for the other changes.

(In reply to Paolo Carlini from comment #20) > Thus I clearly realize something: if we do better for mpfr, the issue > remains that if we do *not* optimize and constants are not propagated, we > issue library calls, which, for command line switches like -fno-math-errno > -funsafe-math-optimizations give incorrect nan result. Where does this nan come from? Do you mean that log(0) from the library gives nan on your machine instead of the correct -inf? For any 0 (as 0 is signed)? Are you sure that there is really a library call?

Thanks. I'm sending to the mailing list an updated version, which only proceeds with do_mpfr_arg1 when !flag_trapping_math && !flag_errno_math. About the -nan, yes by now we know well that if we don't use mpfr at compile time, __builtin_expl(-inf) for -fno-math-errno -funsafe-math-optimizations returns -nan. I think Uros has an explanation why long double is special (in PR57974)

GCC 4.9.0 has been released

GCC 4.9.1 has been released.

GCC 4.9.2 has been released.

GCC 4.9.3 has been released.