[PATCH] Fold sqrt comparisons against constants

[Geoff Keating]

Roger Sayle <roger@www.eyesopen.com> writes:

> This patch performs constant folding of comparisons of the result
> of sqrt, sqrtf, and sqrtl against floating point constants.  This
> optimization can often be applied in scientific software, where
> the distance between two points is compared against a threshold.
> I know of one structural biology application where the single test
> "sqrt(x) < 9.0" is the critical path.  The patch below transforms
> this into the equivalent "(x>=0.0) && (x<81.0)" with -ffast-math.
> 
> The following patch has been tested on i686-pc-linux-gnu with a
> full "make bootstrap", all languages except Ada and treelang, and
> regression tested with a top-level "make -k check" with no new
> failures.
> 
> Ok for mainline?

This is a very good patch, but I think it'd be even better if it was
rearranged a bit.

The cases are as follows:

1. -funsafe-math-optimisations + -ffinite-math-only  (for instance -ffast-math)
   In this case, you can just unconditionally convert 

   sqrt (x) relop y
   to
   x relop y * y

   and perhaps further into 'true' or 'false'.  We don't care about
   overflows or small loss of accuracy with this flag combination.  We
   do, however care that if y*y overflows, then the correct behaviour
   occurs.  (In such a case, the answer will always be 'true' or
   'false', depending on relop and the sign of y).

2. When -fno-finite-math-only is specified
   The important cases to consider are when x < 0, isnan(x), and x ==
   +Inf.  For instance,

   sqrt (x) < DBL_MAX

   is equivalent to 'x >= 0 && x != +Inf'.

   (You also need to consider 'sqrt (x) < nan()', but I think that probably
   Just Works in your previous patch.)

3. When -fno-unsafe-math-optimisations is specified

   It's always possible to remove the sqrt() operation, but it's not
   nearly as simple as just taking 'x relop y * y'.  sqrt() is
   guaranteed to be monotonic, but it compresses its range, so for any
   output value 'y', there will be a range x_0 ... x_1 for which
   sqrt(x) == y.  The problem is that x_0 and x_1 vary depending on
   the rounding mode in use, which GCC could know at compile time but
   at present doesn't.  So I'd suggest skipping this case for now.

> +       else if (code == GT_EXPR || code == GE_EXPR)
> + 	{
> + 	  REAL_VALUE_TYPE c2;
> +
> + 	  /* sqrt(x) > c is the same as x > c*c.  */
> + 	  REAL_ARITHMETIC (c2, MULT_EXPR, c, c);
> + 	  return fold (build (code, type, arg,
> + 			      build_real (TREE_TYPE (arg), c2)));
> + 	}

So this code needs to handle overflow...

> +       else if (code == LT_EXPR || code == LE_EXPR)
> + 	{
> + 	  /* sqrt(x) < c is the same as x >= 0 && x < c*c.  */
> + 	  if ((*lang_hooks.decls.global_bindings_p) () == 0
> + 	      && ! contains_placeholder_p (arg))
> + 	    {
> + 	      REAL_VALUE_TYPE c2;
> +
> + 	      arg = save_expr (arg);
> + 	      REAL_ARITHMETIC (c2, MULT_EXPR, c, c);
> + 	      return fold (build (TRUTH_ANDIF_EXPR, type,
> + 				  fold (build (GE_EXPR, type, arg,
> + 					       build_real (TREE_TYPE (arg),
> + 							   dconst0))),
> + 				  fold (build (code, type, arg,
> + 					       build_real (TREE_TYPE (arg),
> + 							   c2)))));
> + 	    }
> + 	}

... and this code needs to handle overflow, but doesn't need to check
for 'x >= 0' if ! HONOUR_NANS.

-- 
- Geoffrey Keating <geoffk@geoffk.org>