Fix up sqrt(x) < c and sqrt(x) >= c match.pd folding (PR tree-optimization/91734)

Hi!

As mentioned in the PR, the sqrt (x) < c optimization into x < c*c
sometimes breaks the boundary case, if c2=c*c is inexact then in some cases
we need to optimize it into x <= c*c rather than x < c*c.  The original
bugreport is when c is small and c2 is 0.0, then obviously we need <= 0.0
rather than < 0.0, but the testcase includes another example where it makes
a difference, plus has a >= testcase too.

Bootstrapped/regtested on powerpc64le-linux, ok for trunk?

2019-09-13  Jakub Jelinek  <jakub@redhat.com>

	PR tree-optimization/91734
	* generic-match-head.c: Include fold-const-call.h.
	* match.pd (sqrt(x) < c, sqrt(x) >= c): Check the boundary value and
	in case inexact computation of c*c affects comparison of the boundary,
	turn LT_EXPR into LE_EXPR or GE_EXPR into GT_EXPR.

	* gcc.dg/pr91734.c: New test.


	Jakub

On Sat, 14 Sep 2019, Jakub Jelinek wrote:

> Hi!
> 
> As mentioned in the PR, the sqrt (x) < c optimization into x < c*c
> sometimes breaks the boundary case, if c2=c*c is inexact then in some cases
> we need to optimize it into x <= c*c rather than x < c*c.  The original
> bugreport is when c is small and c2 is 0.0, then obviously we need <= 0.0
> rather than < 0.0, but the testcase includes another example where it makes
> a difference, plus has a >= testcase too.
> 
> Bootstrapped/regtested on powerpc64le-linux, ok for trunk?

I was hoping Joseph might chime in here...  anyway, does this assume
round-to-nearest or does it work with round to +-Inf as well?  I
realize this all is under flag_unsafe_math_optimizations, but
this flag is notoriously underspecified...  So the question is
whether we should disable the transform if c*c isn't exact and
flag_rounding_math?  The transform also doesn't seem to guard
against isnan (c) (-funsafe-math-optimizations sets
-fno-trapping-math and -fno-signed-zeros but not -ffinite-math-only
or disables itself on -frounding-math)

Otherwise the patch looks OK to me.

Thanks,
Richard.


> 2019-09-13  Jakub Jelinek  <jakub@redhat.com>
> 
> 	PR tree-optimization/91734
> 	* generic-match-head.c: Include fold-const-call.h.
> 	* match.pd (sqrt(x) < c, sqrt(x) >= c): Check the boundary value and
> 	in case inexact computation of c*c affects comparison of the boundary,
> 	turn LT_EXPR into LE_EXPR or GE_EXPR into GT_EXPR.
> 
> 	* gcc.dg/pr91734.c: New test.
> 
> --- gcc/generic-match-head.c.jj	2019-07-20 21:02:09.296821929 +0200
> +++ gcc/generic-match-head.c	2019-09-12 10:52:33.091366624 +0200
> @@ -29,6 +29,7 @@ along with GCC; see the file COPYING3.
>  #include "cgraph.h"
>  #include "vec-perm-indices.h"
>  #include "fold-const.h"
> +#include "fold-const-call.h"
>  #include "stor-layout.h"
>  #include "tree-dfa.h"
>  #include "builtins.h"
> --- gcc/match.pd.jj	2019-09-11 21:50:54.933504293 +0200
> +++ gcc/match.pd	2019-09-12 11:12:11.150987786 +0200
> @@ -3541,56 +3541,71 @@ (define_operator_list COND_TERNARY
>  	  if x is negative or NaN.  Due to -funsafe-math-optimizations,
>  	  the results for other x follow from natural arithmetic.  */
>         (cmp @0 @1)))
> -     (if (cmp == GT_EXPR || cmp == GE_EXPR)
> +     (if (cmp == LT_EXPR || cmp == LE_EXPR || cmp == GT_EXPR || cmp == GE_EXPR)
>        (with
>         {
> -         REAL_VALUE_TYPE c2;
> +	 REAL_VALUE_TYPE c2;
> +	 enum tree_code ncmp = cmp;
>  	 real_arithmetic (&c2, MULT_EXPR,
>  			  &TREE_REAL_CST (@1), &TREE_REAL_CST (@1));
>  	 real_convert (&c2, TYPE_MODE (TREE_TYPE (@0)), &c2);
> +	 /* See PR91734: if c2 is inexact and sqrt(c2) < c (or sqrt(c2) >= c),
> +	    then change LT_EXPR into LE_EXPR or GE_EXPR into GT_EXPR.  */
> +	 if ((cmp == LT_EXPR || cmp == GE_EXPR) && !REAL_VALUE_ISINF (c2))
> +	   {
> +	     tree c3 = fold_const_call (CFN_SQRT, TREE_TYPE (@0),
> +					build_real (TREE_TYPE (@0), c2));
> +	     if (c3 == NULL_TREE || TREE_CODE (c3) != REAL_CST)
> +	       ncmp = ERROR_MARK;
> +	     else if (real_less (&TREE_REAL_CST (c3), &TREE_REAL_CST (@1)))
> +	       ncmp = cmp == LT_EXPR ? LE_EXPR : GT_EXPR;
> +	   }
>         }
> -       (if (REAL_VALUE_ISINF (c2))
> -	/* sqrt(x) > y is x == +Inf, when y is very large.  */
> -	(if (HONOR_INFINITIES (@0))
> -	 (eq @0 { build_real (TREE_TYPE (@0), c2); })
> -	 { constant_boolean_node (false, type); })
> -	/* sqrt(x) > c is the same as x > c*c.  */
> -	(cmp @0 { build_real (TREE_TYPE (@0), c2); }))))
> -     (if (cmp == LT_EXPR || cmp == LE_EXPR)
> -      (with
> -       {
> -       	 REAL_VALUE_TYPE c2;
> -	 real_arithmetic (&c2, MULT_EXPR,
> -			  &TREE_REAL_CST (@1), &TREE_REAL_CST (@1));
> -	 real_convert (&c2, TYPE_MODE (TREE_TYPE (@0)), &c2);
> -       }
> -       (if (REAL_VALUE_ISINF (c2))
> -        (switch
> -	 /* sqrt(x) < y is always true, when y is a very large
> -	    value and we don't care about NaNs or Infinities.  */
> -	 (if (! HONOR_NANS (@0) && ! HONOR_INFINITIES (@0))
> -	  { constant_boolean_node (true, type); })
> -	 /* sqrt(x) < y is x != +Inf when y is very large and we
> -	    don't care about NaNs.  */
> -	 (if (! HONOR_NANS (@0))
> -	  (ne @0 { build_real (TREE_TYPE (@0), c2); }))
> -	 /* sqrt(x) < y is x >= 0 when y is very large and we
> -	    don't care about Infinities.  */
> -	 (if (! HONOR_INFINITIES (@0))
> -	  (ge @0 { build_real (TREE_TYPE (@0), dconst0); }))
> -	 /* sqrt(x) < y is x >= 0 && x != +Inf, when y is large.  */
> -	 (if (GENERIC)
> -	  (truth_andif
> -	   (ge @0 { build_real (TREE_TYPE (@0), dconst0); })
> -	   (ne @0 { build_real (TREE_TYPE (@0), c2); }))))
> -	/* sqrt(x) < c is the same as x < c*c, if we ignore NaNs.  */
> -	(if (! HONOR_NANS (@0))
> -	 (cmp @0 { build_real (TREE_TYPE (@0), c2); })
> -	 /* sqrt(x) < c is the same as x >= 0 && x < c*c.  */
> -	 (if (GENERIC)
> -	  (truth_andif
> -	   (ge @0 { build_real (TREE_TYPE (@0), dconst0); })
> -	   (cmp @0 { build_real (TREE_TYPE (@0), c2); })))))))))
> +       (if (cmp == GT_EXPR || cmp == GE_EXPR)
> +	(if (REAL_VALUE_ISINF (c2))
> +	 /* sqrt(x) > y is x == +Inf, when y is very large.  */
> +	 (if (HONOR_INFINITIES (@0))
> +	  (eq @0 { build_real (TREE_TYPE (@0), c2); })
> +	  { constant_boolean_node (false, type); })
> +	 /* sqrt(x) > c is the same as x > c*c.  */
> +	 (if (ncmp != ERROR_MARK)
> +	  (if (ncmp == GE_EXPR)
> +	   (ge @0 { build_real (TREE_TYPE (@0), c2); })
> +	   (gt @0 { build_real (TREE_TYPE (@0), c2); }))))
> +	/* else if (cmp == LT_EXPR || cmp == LE_EXPR)  */
> +	(if (REAL_VALUE_ISINF (c2))
> +	 (switch
> +	  /* sqrt(x) < y is always true, when y is a very large
> +	     value and we don't care about NaNs or Infinities.  */
> +	  (if (! HONOR_NANS (@0) && ! HONOR_INFINITIES (@0))
> +	   { constant_boolean_node (true, type); })
> +	  /* sqrt(x) < y is x != +Inf when y is very large and we
> +	     don't care about NaNs.  */
> +	  (if (! HONOR_NANS (@0))
> +	   (ne @0 { build_real (TREE_TYPE (@0), c2); }))
> +	  /* sqrt(x) < y is x >= 0 when y is very large and we
> +	     don't care about Infinities.  */
> +	  (if (! HONOR_INFINITIES (@0))
> +	   (ge @0 { build_real (TREE_TYPE (@0), dconst0); }))
> +	  /* sqrt(x) < y is x >= 0 && x != +Inf, when y is large.  */
> +	  (if (GENERIC)
> +	   (truth_andif
> +	    (ge @0 { build_real (TREE_TYPE (@0), dconst0); })
> +	    (ne @0 { build_real (TREE_TYPE (@0), c2); }))))
> +	 /* sqrt(x) < c is the same as x < c*c, if we ignore NaNs.  */
> +	 (if (ncmp != ERROR_MARK && ! HONOR_NANS (@0))
> +	  (if (ncmp == LT_EXPR)
> +	   (lt @0 { build_real (TREE_TYPE (@0), c2); })
> +	   (le @0 { build_real (TREE_TYPE (@0), c2); }))
> +	  /* sqrt(x) < c is the same as x >= 0 && x < c*c.  */
> +	  (if (ncmp != ERROR_MARK && GENERIC)
> +	   (if (ncmp == LT_EXPR)
> +	    (truth_andif
> +	     (ge @0 { build_real (TREE_TYPE (@0), dconst0); })
> +	     (lt @0 { build_real (TREE_TYPE (@0), c2); }))
> +	    (truth_andif
> +	     (ge @0 { build_real (TREE_TYPE (@0), dconst0); })
> +	     (le @0 { build_real (TREE_TYPE (@0), c2); })))))))))))
>     /* Transform sqrt(x) cmp sqrt(y) -> x cmp y.  */
>     (simplify
>      (cmp (sq @0) (sq @1))
> --- gcc/testsuite/gcc.dg/pr91734.c.jj	2019-09-12 10:52:33.094366596 +0200
> +++ gcc/testsuite/gcc.dg/pr91734.c	2019-09-12 10:49:10.000000000 +0200
> @@ -0,0 +1,34 @@
> +/* PR tree-optimization/91734 */
> +/* { dg-do run } */
> +/* { dg-add-options ieee } */
> +/* { dg-additional-options "-ffast-math -O2 -std=gnu99" } */
> +
> +__attribute__((noipa)) int
> +foo (float x)
> +{
> +  return __builtin_sqrtf (x) < __FLT_MIN__;
> +}
> +
> +__attribute__((noipa)) int
> +bar (float x)
> +{
> +  return __builtin_sqrtf (x) < 0x1.2dd3d0p-65f;
> +}
> +
> +__attribute__((noipa)) int
> +baz (float x)
> +{
> +  return __builtin_sqrtf (x) >= 0x1.2dd3d0p-65f;
> +}
> +
> +int
> +main ()
> +{
> +  if (!foo (0.0f))
> +    __builtin_abort ();
> +  if (!bar (0x1.63dbc0p-130f))
> +    __builtin_abort ();
> +  if (baz (0x1.63dbc0p-130f))
> +    __builtin_abort ();
> +  return 0;
> +}
> 
> 	Jakub
>

On Sat, 14 Sep 2019, Jakub Jelinek wrote:

> Hi!
> 
> As mentioned in the PR, the sqrt (x) < c optimization into x < c*c
> sometimes breaks the boundary case, if c2=c*c is inexact then in some cases
> we need to optimize it into x <= c*c rather than x < c*c.  The original

And in some cases the boundary case is wrong in the other direction and we 
need to optimize it into x < nextdown (c*c).  For example, say c == 
0x1.001002p+0f; (double) c * c == 0x1.002005004004p+0, and rounding to 
nearest float gives c * c == 0x1.002006p+0f, so nextdownf (c * c) == 
0x1.002004p+0f; in double precision, sqrt (0x1.002004p+0f) == 
0x1.0010017fe8006p+0, and rounding to nearest float gives c again.  So x < 
c*c does not in fact imply sqrt (x) < c in that case.