This is the mail archive of the gcc-patches@gcc.gnu.org mailing list for the GCC project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[PATCH] Evaluate pow(x,n) at compile-time (take 2)


The following patch is a revision of my previous patch to evaluate
pow(x,n), for integer constant values of n, at compile time.
http://gcc.gnu.org/ml/gcc-patches/2003-04/msg00627.html

I believe this revision addresses most concerns raised over the
original patch, primarily by Geoff Keating.  To this end, I've
since modified real.c to allow the precision of results to be
tracked http://gcc.gnu.org/ml/gcc-patches/2003-04/msg01543.html
This now allows us to know when "pow(x,n)" is guaranteed to be
exact, 0 ulp error, and when it may have been rounded.  This is
overly pessimistic as correctly rounded results are still perfect,
i.e. less than 0.5 ulp.

With this infrastructure real_powi now returns a flag indicating
whether the result is exact or not, in which case we can safely
apply the transformation even without "unsafe" math optimizations.

I still argue that with -funsafe-math-optimizations we allow
pow(0.1,3) to be equivalent to 0.1*0.1*0.1 as is performed by
other optimizing compilers, and even the g77 front-end.  Without
this gcc and g++ will be slower than FORTRAN on scientific and
numerical codes that make heavy use of exponentiation.  After
all this is exactly what "unsafe math optimizations" means.

I'm happy to post my back of the envelope error analysis of the
implementation should anyone be interested.  My numerical analysis
may be rusty but I believe that even with unsafe math optimizations
we're always much less than 1 ulp for IEEE double, as required
for the system library.


The following patch has been tested on i686-pc-linux-gnu with a
complete "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?


2003-05-04  Roger Sayle  <roger@eyesopen.com>

	* real.c (real_powi): New function to calculate the value of
	a real raised to an integer power, i.e. pow(x,n) for int n.
	(real_sqrt): Convert to using the faster do_add, do_multiply
	and do_divide API for consistency with the rest of real.c.
	* real.h (real_powi): Prototype here.
	* builtins.c (fold_builtin):  Avoid local variable mode when
	evaluating sqrt at compile time.  Attempt to evaluate pow at
	compile-time, by checking for an integral exponent.

	* gcc.dg/builtins-14.c: New test case.


Index: real.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/real.c,v
retrieving revision 1.117
diff -c -3 -p -r1.117 real.c
*** real.c	23 Apr 2003 02:46:03 -0000	1.117
--- real.c	4 May 2003 02:24:26 -0000
*************** real_sqrt (r, mode, x)
*** 4606,4612 ****

    if (!init)
      {
!       real_arithmetic (&halfthree, PLUS_EXPR, &dconst1, &dconsthalf);
        init = true;
      }

--- 4606,4612 ----

    if (!init)
      {
!       do_add (&halfthree, &dconst1, &dconsthalf, 0);
        init = true;
      }

*************** real_sqrt (r, mode, x)
*** 4618,4628 ****
    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, &dconsthalf);
!       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))
--- 4618,4628 ----
    for (iter = 0; iter < 16; iter++)
      {
        /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x).  */
!       do_multiply (&t, x, &i);
!       do_multiply (&h, &t, &i);
!       do_multiply (&t, &h, &dconsthalf);
!       do_add (&h, &halfthree, &t, 1);
!       do_multiply (&t, &i, &h);

        /* Check for early convergence.  */
        if (iter >= 6 && real_identical (&i, &t))
*************** real_sqrt (r, mode, x)
*** 4633,4648 ****
      }

    /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)).  */
!   real_arithmetic (&t, MULT_EXPR, x, &i);
!   real_arithmetic (&h, MULT_EXPR, &t, &i);
!   real_arithmetic (&i, MINUS_EXPR, &dconst1, &h);
!   real_arithmetic (&h, MULT_EXPR, &t, &i);
!   real_arithmetic (&i, MULT_EXPR, &dconsthalf, &h);
!   real_arithmetic (&h, PLUS_EXPR, &t, &i);

    /* ??? We need a Tuckerman test to get the last bit.  */

    real_convert (r, mode, &h);
    return true;
  }

--- 4633,4704 ----
      }

    /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)).  */
!   do_multiply (&t, x, &i);
!   do_multiply (&h, &t, &i);
!   do_add (&i, &dconst1, &h, 1);
!   do_multiply (&h, &t, &i);
!   do_multiply (&i, &dconsthalf, &h);
!   do_add (&h, &t, &i, 0);

    /* ??? We need a Tuckerman test to get the last bit.  */

    real_convert (r, mode, &h);
    return true;
+ }
+
+ /* Calculate X raised to the integer exponent N in mode MODE and store
+    the result in R.  Return true if the result may be inexact due to
+    loss of precision.  The algorithm is the classic "left-to-right binary
+    method" described in section 4.6.3 of Donald Knuth's "Seminumerical
+    Algorithms", "The Art of Computer Programming", Volume 2.  */
+
+ bool
+ real_powi (r, mode, x, n)
+      REAL_VALUE_TYPE *r;
+      enum machine_mode mode;
+      const REAL_VALUE_TYPE *x;
+      HOST_WIDE_INT n;
+ {
+   unsigned HOST_WIDE_INT bit;
+   REAL_VALUE_TYPE t;
+   bool inexact = false;
+   bool init = false;
+   bool neg;
+   int i;
+
+   if (n == 0)
+     {
+       *r = dconst1;
+       return false;
+     }
+   else if (n < 0)
+     {
+       /* Don't worry about overflow, from now on n is unsigned.  */
+       neg = true;
+       n = -n;
+     }
+   else
+     neg = false;
+
+   t = *x;
+   bit = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
+   for (i = 0; i < HOST_BITS_PER_WIDE_INT; i++)
+     {
+       if (init)
+ 	{
+ 	  inexact |= do_multiply (&t, &t, &t);
+ 	  if (n & bit)
+ 	    inexact |= do_multiply (&t, &t, x);
+ 	}
+       else if (n & bit)
+ 	init = true;
+       bit >>= 1;
+     }
+
+   if (neg)
+     inexact |= do_divide (&t, &dconst1, &t);
+
+   real_convert (r, mode, &t);
+   return inexact;
  }

Index: real.h
===================================================================
RCS file: /cvs/gcc/gcc/gcc/real.h,v
retrieving revision 1.64
diff -c -3 -p -r1.64 real.h
*** real.h	1 Apr 2003 21:45:27 -0000	1.64
--- real.h	4 May 2003 02:24:26 -0000
*************** extern bool real_sqrt			PARAMS ((REAL_VA
*** 365,368 ****
--- 365,374 ----
  						 enum machine_mode,
  						 const REAL_VALUE_TYPE *));

+ /* Calculate R as X raised to the integer exponent N in mode MODE. */
+ extern bool real_powi			PARAMS ((REAL_VALUE_TYPE *,
+ 						 enum machine_mode,
+ 						 const REAL_VALUE_TYPE *,
+ 						 HOST_WIDE_INT));
+
  #endif /* ! GCC_REAL_H */
Index: builtins.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/builtins.c,v
retrieving revision 1.194
diff -c -3 -p -r1.194 builtins.c
*** builtins.c	3 May 2003 14:25:20 -0000	1.194
--- builtins.c	4 May 2003 02:24:27 -0000
*************** fold_builtin (exp)
*** 4967,4978 ****
  	  if (TREE_CODE (arg) == REAL_CST
  	      && ! TREE_CONSTANT_OVERFLOW (arg))
  	    {
- 	      enum machine_mode mode;
  	      REAL_VALUE_TYPE r, x;

  	      x = TREE_REAL_CST (arg);
! 	      mode = TYPE_MODE (type);
! 	      if (real_sqrt (&r, mode, &x)
  		  || (!flag_trapping_math && !flag_errno_math))
  		return build_real (type, r);
  	    }
--- 4967,4976 ----
  	  if (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 (type), &x)
  		  || (!flag_trapping_math && !flag_errno_math))
  		return build_real (type, r);
  	    }
*************** fold_builtin (exp)
*** 5183,5188 ****
--- 5181,5208 ----
  		    {
  		      tree arglist = build_tree_list (NULL_TREE, arg0);
  		      return build_function_call_expr (sqrtfn, arglist);
+ 		    }
+ 		}
+
+ 	      /* Attempt to evaluate pow at compile-time.  */
+ 	      if (TREE_CODE (arg0) == REAL_CST
+ 		  && ! TREE_CONSTANT_OVERFLOW (arg0))
+ 		{
+ 		  REAL_VALUE_TYPE cint;
+ 		  HOST_WIDE_INT n;
+
+ 		  n = real_to_integer(&c);
+ 		  real_from_integer (&cint, VOIDmode, n,
+ 				     n < 0 ? -1 : 0, 0);
+ 		  if (real_identical (&c, &cint))
+ 		    {
+ 		      REAL_VALUE_TYPE x;
+ 		      bool inexact;
+
+ 		      x = TREE_REAL_CST (arg0);
+ 		      inexact = real_powi (&x, TYPE_MODE (type), &x, n);
+ 		      if (flag_unsafe_math_optimizations || !inexact)
+ 			return build_real (type, x);
  		    }
  		}
  	    }


/* Copyright (C) 2003 Free Software Foundation.

   Check that constant folding of built-in math functions doesn't
   break anything and produces the expected results.

   Written by Roger Sayle, 9th April 2003.  */

/* { dg-do link } */
/* { dg-options "-O2" } */

extern void link_error(void);

extern double pow(double,double);


int main()
{
  if (pow (2.0, 3.0) != 8.0)
    link_error ();

  if (pow (2.0, -3.0) != 0.125)
    link_error ();

  return 0;
}

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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]