From e82a312b55f88fc3d035e0c24409c71e043e9633 Mon Sep 17 00:00:00 2001
From: Roger Sayle <roger@eyesopen.com>
Date: Tue, 6 May 2003 03:14:10 +0000
Subject: real.c (real_powi): New function to calculate the value of a real
 raised to an integer power, i.e.

	* 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.

From-SVN: r66515
---
 gcc/real.c | 80 ++++++++++++++++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 68 insertions(+), 12 deletions(-)

(limited to 'gcc/real.c')

diff --git a/gcc/real.c b/gcc/real.c
index db658fd..16cd02a 100644
--- a/gcc/real.c
+++ b/gcc/real.c
@@ -4606,7 +4606,7 @@ real_sqrt (r, mode, x)
 
   if (!init)
     {
-      real_arithmetic (&halfthree, PLUS_EXPR, &dconst1, &dconsthalf);
+      do_add (&halfthree, &dconst1, &dconsthalf, 0);
       init = true;
     }
 
@@ -4618,11 +4618,11 @@ real_sqrt (r, mode, x)
   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);
+      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))
@@ -4633,12 +4633,12 @@ real_sqrt (r, mode, x)
     }
 
   /* 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);
+  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.  */
 
@@ -4646,3 +4646,59 @@ real_sqrt (r, mode, x)
   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;
+}
+
-- 
cgit v1.1