[PATCH] Constant fold sqrt at compile time

[Roger Sayle <roger at eyesopen dot com>]

I was saving this patch until after 3.3 branched, but with the
recent interest in the "paranoia" thread on the gcc list, I thought
I'd post it for review whilst we have the attention of the number
theorists.

This patch performs constant folding and RTL simplification of
sqrt, sqrtf and sqrtl with a constant argument at compile time.

I'd like to thank Brad Lucier for his helpful comments.  It was
his suggestion that real_sqrt take the desired mode of the result
such that an implementation can make use of it to avoid double
rouding and achieve 0.5 ulp.  The current implementation however
doesn't make use of this, and computes a 160-bit square root to
about 2 ulps, which is then truncated to the precision of the
result.  For this reason, evaluation at compile time is currently
guarded by -funsafe-math-optimizations/-ffast-math.

My goal is that this patch might serve as an example as to how
to add similar functionality for exp, log, sin and cos (and
perhaps pow) to GCC.  Improved implementations of sqrt are also
welcome.  Brad hinted that if accepted he may be able to post
a solution to achieve 0.5 ulp, in which case we could enable
this optimization permanently.

Finally, although I've placed the implementation of real_sqrt
in real.c, it only uses the external real.h API, allowing it
to be moved to a separate file if required.

The following patch against the gcc-3_4-basic-improvements-branch
has been tested by a full "make bootstrap" and "make -k check"
all languages except Ada and treelang, on i686-pc-linux-gnu with
no new regressions.  It is triggered in the testsuite by my earlier
tests to check that "sqrt(0.0) == 0.0" and "sqrt(1.0) == 1.0",
and by hand with several random values.

Ok for gcc-3_4-basic-improvements-branch?


2002-10-16  Roger Sayle  <roger@eyesopen.com>

	* real.c (real_sqrt): New function to calculate square roots.
	* real.h (real_sqrt): Add function prototype.
	* builtins.c (fold_builtin): Fold sqrt of constant argument.
	* simplify-rtx.c (simplify_unary_operation): Simplify sqrt
	of constant argument.


Index: real.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/real.c,v
retrieving revision 1.75.4.5
diff -c -3 -p -r1.75.4.5 real.c
*** real.c	15 Oct 2002 01:32:58 -0000	1.75.4.5
--- real.c	17 Oct 2002 05:25:24 -0000
*************** const struct real_format *real_format_fo
*** 4106,4108 ****
--- 4106,4182 ----
    NULL,				/* XFmode */
    &ieee_quad_format		/* TFmode */
  };
+
+
+ /* Calculate the square root of X in mode MODE, and store the result
+    in R.  Returns true if this doesn't trap and doesn't set errno.  */
+
+ bool
+ real_sqrt (r, mode, x)
+      REAL_VALUE_TYPE *r;
+      enum machine_mode mode;
+      const REAL_VALUE_TYPE *x;
+ {
+   static REAL_VALUE_TYPE halfthree;
+   static REAL_VALUE_TYPE half;
+   static bool init = false;
+   REAL_VALUE_TYPE h, t, i;
+   int iter, exp;
+
+   /* sqrt(-0.0) is -0.0.  */
+   if (real_isnegzero (x))
+     {
+       *r = *x;
+       return false;
+     }
+
+   /* Negative arguments return NaN.  */
+   if (real_isneg (x))
+     {
+       /* Mode is ignored for canonical NaN.  */
+       real_nan (r, "", 1, SFmode);
+       return false;
+     }
+
+   /* Infinity and NaN return themselves.  */
+   if (real_isinf (x) || real_isnan (x))
+     {
+       *r = *x;
+       return false;
+     }
+
+   if (!init)
+     {
+       real_arithmetic (&half, RDIV_EXPR, &dconst1, &dconst2);
+       real_arithmetic (&halfthree, PLUS_EXPR, &dconst1, &half);
+       init = true;
+     }
+
+   /* Initial guess for reciprocal sqrt, i.  */
+   exp = real_exponent (x);
+   real_ldexp (&i, &dconst1, -exp/2);
+
+   /* Newton's iteration for reciprocal sqrt, i */
+   for (iter = 0; iter < 16; iter++)
+     {
+       /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x).  */
+       real_arithmetic (&t, MULT_EXPR, x, &i);
+       real_arithmetic (&h, MULT_EXPR, &t, &i);
+       real_arithmetic (&t, MULT_EXPR, &h, &half);
+       real_arithmetic (&h, MINUS_EXPR, &halfthree, &t);
+       real_arithmetic (&t, MULT_EXPR, &i, &h);
+
+       /* Check for early convergence.  */
+       if (iter >= 6 && real_identical (&i, &t))
+ 	break;
+
+       /* ??? Unroll loop to avoid copying.  */
+       i = t;
+     }
+
+   /* Final result is original value, x, multiplied by i.  */
+   real_arithmetic (&t, MULT_EXPR, &i, x);
+   real_convert (r, mode, &t);
+   return true;
+ }
+
Index: real.h
===================================================================
RCS file: /cvs/gcc/gcc/gcc/real.h,v
retrieving revision 1.44.8.4
diff -c -3 -p -r1.44.8.4 real.h
*** real.h	15 Oct 2002 01:32:59 -0000	1.44.8.4
--- real.h	17 Oct 2002 05:25:24 -0000
*************** extern bool exact_real_inverse	PARAMS ((
*** 348,352 ****
--- 348,356 ----
  /* In tree.c: wrap up a REAL_VALUE_TYPE in a tree node.  */
  extern tree build_real			PARAMS ((tree, REAL_VALUE_TYPE));

+ /* Calculate R as the square root of X in the given machine mode.  */
+ extern bool real_sqrt			PARAMS ((REAL_VALUE_TYPE *,
+ 						 enum machine_mode,
+ 						 const REAL_VALUE_TYPE *));

  #endif /* ! GCC_REAL_H */
Index: builtins.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/builtins.c,v
retrieving revision 1.159.4.3
diff -c -3 -p -r1.159.4.3 builtins.c
*** builtins.c	17 Sep 2002 22:58:36 -0000	1.159.4.3
--- builtins.c	17 Oct 2002 05:25:25 -0000
*************** fold_builtin (exp)
*** 4241,4249 ****
  	  enum built_in_function fcode;
  	  tree arg = TREE_VALUE (arglist);

! 	  /* Optimize sqrt(0.0) = 0.0 and sqrt(1.0) = 1.0.  */
! 	  if (real_zerop (arg) || real_onep (arg))
! 	    return arg;

  	  /* Optimize sqrt(exp(x)) = exp(x/2.0).  */
  	  fcode = builtin_mathfn_code (arg);
--- 4241,4258 ----
  	  enum built_in_function fcode;
  	  tree arg = TREE_VALUE (arglist);

! 	  /* Optimize sqrt of constant value.  */
! 	  if (flag_unsafe_math_optimizations
! 	      && TREE_CODE (arg) == REAL_CST
! 	      && ! TREE_CONSTANT_OVERFLOW (arg))
! 	    {
! 	      REAL_VALUE_TYPE r, x;
!
! 	      x = TREE_REAL_CST (arg);
! 	      if (real_sqrt (&r, TYPE_MODE (TREE_TYPE (arg)), &x)
! 		  || (!flag_trapping_math && !flag_errno_math))
! 		return build_real (TREE_TYPE (arg), r);
! 	    }

  	  /* Optimize sqrt(exp(x)) = exp(x/2.0).  */
  	  fcode = builtin_mathfn_code (arg);
Index: simplify-rtx.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/simplify-rtx.c,v
retrieving revision 1.118.4.5
diff -c -3 -p -r1.118.4.5 simplify-rtx.c
*** simplify-rtx.c	15 Oct 2002 01:33:01 -0000	1.118.4.5
--- simplify-rtx.c	17 Oct 2002 05:25:25 -0000
*************** simplify_unary_operation (code, mode, op
*** 558,572 ****
    else if (GET_CODE (trueop) == CONST_DOUBLE
  	   && GET_MODE_CLASS (mode) == MODE_FLOAT)
      {
!       REAL_VALUE_TYPE d;
        REAL_VALUE_FROM_CONST_DOUBLE (d, trueop);

        switch (code)
  	{
  	case SQRT:
! 	  /* We don't attempt to optimize this.  */
! 	  return 0;
!
  	case ABS:
  	  d = REAL_VALUE_ABS (d);
  	  break;
--- 558,576 ----
    else if (GET_CODE (trueop) == CONST_DOUBLE
  	   && GET_MODE_CLASS (mode) == MODE_FLOAT)
      {
!       REAL_VALUE_TYPE d, t;
        REAL_VALUE_FROM_CONST_DOUBLE (d, trueop);

        switch (code)
  	{
  	case SQRT:
! 	  if (! flag_unsafe_math_optimizations)
! 	    return 0;
! 	  if (! real_sqrt (&t, mode, &d)
! 	      && (flag_trapping_math || flag_errno_math))
! 	    return 0;
! 	  d = t;
! 	  break;
  	case ABS:
  	  d = REAL_VALUE_ABS (d);
  	  break;

Roger
--
Roger Sayle,                         E-mail: roger@eyesopen.com
OpenEye Scientific Software,         WWW: http://www.eyesopen.com/
Suite 1107, 3600 Cerrillos Road,     Tel: (+1) 505-473-7385
Santa Fe, New Mexico, 87507.         Fax: (+1) 505-473-0833