Bug 57994 - Constant folding of infinity
Summary: Constant folding of infinity
Status: ASSIGNED
Alias: None
Product: gcc
Classification: Unclassified
Component: tree-optimization (show other bugs)
Version: 4.9.0
: P3 normal
Target Milestone: ---
Assignee: Paolo Carlini
URL:
Keywords: missed-optimization
Depends on:
Blocks:
 
Reported: 2013-07-26 15:13 UTC by Marc Glisse
Modified: 2015-06-26 20:49 UTC (History)
3 users (show)

See Also:
Host:
Target:
Build:
Known to work:
Known to fail:
Last reconfirmed: 2013-07-27 00:00:00


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Marc Glisse 2013-07-26 15:13:47 UTC
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).
Comment 1 Paolo Carlini 2013-07-26 15:16:08 UTC
Thanks Marc.
Comment 2 joseph@codesourcery.com 2013-07-26 21:53:52 UTC
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.
Comment 3 Paolo Carlini 2013-07-27 00:15:14 UTC
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?
Comment 4 Marc Glisse 2013-07-27 07:55:32 UTC
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.
Comment 5 Paolo Carlini 2013-07-27 19:53:49 UTC
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?
Comment 6 Marc Glisse 2013-07-27 20:10:03 UTC
(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.
Comment 7 joseph@codesourcery.com 2013-07-27 20:10:51 UTC
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.
Comment 8 Paolo Carlini 2013-07-27 20:26:47 UTC
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.
Comment 9 Marc Glisse 2013-07-27 20:31:25 UTC
(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.
Comment 10 Paolo Carlini 2013-07-27 20:46:00 UTC
About testing, it would be just matter of extending/updating what Kaveh Ghazi set up when mpfr/mpc came in.
Comment 11 joseph@codesourcery.com 2013-07-28 16:37:55 UTC
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).
Comment 12 Vincent Lefèvre 2013-07-28 17:40:42 UTC
(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).
Comment 13 Vincent Lefèvre 2013-08-01 20:55:18 UTC
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.
Comment 14 Paolo Carlini 2013-08-02 23:23:36 UTC
Let's see.
Comment 15 Vincent Lefèvre 2013-08-27 17:34:41 UTC
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
Comment 16 Paolo Carlini 2013-10-24 12:26:59 UTC
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.
Comment 17 Paolo Carlini 2013-10-24 12:35:53 UTC
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).
Comment 18 Paolo Carlini 2013-10-24 12:40:10 UTC
Now the question is: why we use a library call for log(0) instead of mpfr?
Comment 19 Paolo Carlini 2013-10-24 13:01:38 UTC
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)))
         {
Comment 20 Paolo Carlini 2013-10-24 13:08:02 UTC
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!
Comment 21 Vincent Lefèvre 2013-10-24 15:29:50 UTC
(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.
Comment 22 Vincent Lefèvre 2013-10-24 15:42:52 UTC
(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?
Comment 23 Paolo Carlini 2013-10-24 15:51:21 UTC
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)
Comment 24 Jakub Jelinek 2014-04-22 11:36:30 UTC
GCC 4.9.0 has been released
Comment 25 Jakub Jelinek 2014-07-16 13:29:16 UTC
GCC 4.9.1 has been released.
Comment 26 Jakub Jelinek 2014-10-30 10:39:22 UTC
GCC 4.9.2 has been released.
Comment 27 Jakub Jelinek 2015-06-26 19:57:07 UTC
GCC 4.9.3 has been released.