[PATCH]: PR29335 evaluate transcendentals at compile-time using MPFR [take 3]

["Kaveh R. GHAZI" <ghazi at caip dot rutgers dot edu>]

This is is my third attempt to have GCC evaluate transcendentals at
compile-time using MPFR.  I've tested the patch using the
transcendental accuracy benchmark pointed out by Roger here:
http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00760.html

With some offline help from Roger (and a bit of duct tape) I was able to
get the thing testing GCC manually.  (Getting this test into GCC itself
would involve a lot more work, plus the assignment of the author which I
don't have.)  I'm happy to note that with this revision of my patch, I'm
getting "perfect" results, i.e. 0.5 or less ulp error for all the
functions I've converted so far (sin, cos, tan).

Changes from the last revision are:

1.  Add check for ! TREE_CONSTANT_OVERFLOW (arg).

2.  Ensure that (arg) is not nan or inf.

3.  Use a hex (not decimal) intermediate string in mpfr_from_real().

4.  Use real_convert() to convert the result into the requested mode.

5.  Run the transcendental function in double precision, then round to
    the requested precision and check for exact results.

Let me explain #5.  The previous version of the patch set the MPFR
precision to exactly that of the supplied type.  I.e. for double, it
would use 53.  I did this so that I could see if the result was exact,
then I could fold even if the user requested -frounding-math.  See:
http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00374.html
http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00653.html

When using this precision, for the sin tests I was getting three of
the results with errors of 0.50000 ... 000N, i.e. ever so slightly
above 0.5.  To put that in context, sin has 1257 tests, cos has 1418,
tan has 1596.  So 1254 of the sin tests were perfect, and the cos &
tan results were all 100% perfect.

I started increasing the MPFR precision and the errors got smaller. At
precisio=75, they were all gone.  I didn't see anything magical about 75,
so I picked "supplied precision * 2" as the precision within which mpfr
should do it's transcendental calculation and after that function call I
round to the requested precision.  I check that both operations were exact
before proceeding with -frounding-math.  This seems to do the trick.

On solaris10, the libc sin's worst error was 0.57, cos' was 0.56 and
tan's was 0.51.  That's not bad, but now with my patch these "hard
cases" from the transcendental testsuite all get 0.5 ulp error or less
for constant arguments.

Bootstrapped on sparc-sun-solaris2.10, no regressions.  This patch
relies on the GMP/MPFR configury patch previously posted here:
http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00720.html

Is this one okay for mainline stage1?

		Thanks,
		--Kaveh


2006-10-19  Kaveh R. Ghazi  <ghazi@caip.rutgers.edu>

	PR middle-end/29335
	* builtins.c (fold_builtin_sin, fold_builtin_cos,
	fold_builtin_tan): Fold all constant arguments.  Take a "type"
	argument as necessary.
	(do_mpfr_arg1): New.
	* real.c, real.h (real_from_mpfr, mpfr_from_real): New.

diff -rup orig/egcc-SVN20061016/gcc/builtins.c egcc-SVN20061016/gcc/builtins.c
--- orig/egcc-SVN20061016/gcc/builtins.c	2006-10-15 20:01:40.000000000 -0400
+++ egcc-SVN20061016/gcc/builtins.c	2006-10-19 20:49:27.781859356 -0400
@@ -149,9 +149,9 @@ static tree fold_builtin_sqrt (tree, tre
 static tree fold_builtin_cbrt (tree, tree);
 static tree fold_builtin_pow (tree, tree, tree);
 static tree fold_builtin_powi (tree, tree, tree);
-static tree fold_builtin_sin (tree);
+static tree fold_builtin_sin (tree, tree);
 static tree fold_builtin_cos (tree, tree, tree);
-static tree fold_builtin_tan (tree);
+static tree fold_builtin_tan (tree, tree);
 static tree fold_builtin_atan (tree, tree);
 static tree fold_builtin_trunc (tree, tree);
 static tree fold_builtin_floor (tree, tree);
@@ -204,6 +204,7 @@ static unsigned HOST_WIDE_INT target_s;
 static char target_percent_c[3];
 static char target_percent_s[3];
 static char target_percent_s_newline[4];
+static tree do_mpfr_arg1 (tree, tree, int (*)(mpfr_ptr, mpfr_srcptr, mp_rnd_t));

 /* Return true if NODE should be considered for inline expansion regardless
    of the optimization level.  This means whenever a function is invoked with
@@ -7157,17 +7158,17 @@ fold_builtin_cbrt (tree arglist, tree ty
 /* Fold function call to builtin sin, sinf, or sinl.  Return
    NULL_TREE if no simplification can be made.  */
 static tree
-fold_builtin_sin (tree arglist)
+fold_builtin_sin (tree arglist, tree type)
 {
-  tree arg = TREE_VALUE (arglist);
+  tree arg = TREE_VALUE (arglist), res;

   if (!validate_arglist (arglist, REAL_TYPE, VOID_TYPE))
     return NULL_TREE;

-  /* Optimize sin (0.0) = 0.0.  */
-  if (real_zerop (arg))
-    return arg;
-
+  /* Calculate the result when the argument is a constant.  */
+  if ((res = do_mpfr_arg1 (arg, type, mpfr_sin)))
+    return res;
+
   return NULL_TREE;
 }

@@ -7176,15 +7177,15 @@ fold_builtin_sin (tree arglist)
 static tree
 fold_builtin_cos (tree arglist, tree type, tree fndecl)
 {
-  tree arg = TREE_VALUE (arglist);
+  tree arg = TREE_VALUE (arglist), res;

   if (!validate_arglist (arglist, REAL_TYPE, VOID_TYPE))
     return NULL_TREE;

-  /* Optimize cos (0.0) = 1.0.  */
-  if (real_zerop (arg))
-    return build_real (type, dconst1);
-
+  /* Calculate the result when the argument is a constant.  */
+  if ((res = do_mpfr_arg1 (arg, type, mpfr_cos)))
+    return res;
+
   /* Optimize cos(-x) into cos (x).  */
   if (TREE_CODE (arg) == NEGATE_EXPR)
     {
@@ -7199,18 +7200,18 @@ fold_builtin_cos (tree arglist, tree typ
 /* Fold function call to builtin tan, tanf, or tanl.  Return
    NULL_TREE if no simplification can be made.  */
 static tree
-fold_builtin_tan (tree arglist)
+fold_builtin_tan (tree arglist, tree type)
 {
   enum built_in_function fcode;
-  tree arg = TREE_VALUE (arglist);
+  tree arg = TREE_VALUE (arglist), res;

   if (!validate_arglist (arglist, REAL_TYPE, VOID_TYPE))
     return NULL_TREE;

-  /* Optimize tan(0.0) = 0.0.  */
-  if (real_zerop (arg))
-    return arg;
-
+  /* Calculate the result when the argument is a constant.  */
+  if ((res = do_mpfr_arg1 (arg, type, mpfr_tan)))
+    return res;
+
   /* Optimize tan(atan(x)) = x.  */
   fcode = builtin_mathfn_code (arg);
   if (flag_unsafe_math_optimizations
@@ -8952,7 +8953,7 @@ fold_builtin_1 (tree fndecl, tree arglis
       return fold_builtin_cbrt (arglist, type);

     CASE_FLT_FN (BUILT_IN_SIN):
-      return fold_builtin_sin (arglist);
+      return fold_builtin_sin (arglist, type);

     CASE_FLT_FN (BUILT_IN_COS):
       return fold_builtin_cos (arglist, type, fndecl);
@@ -8977,7 +8978,7 @@ fold_builtin_1 (tree fndecl, tree arglis
       return fold_builtin_logarithm (fndecl, arglist, &dconst10);

     CASE_FLT_FN (BUILT_IN_TAN):
-      return fold_builtin_tan (arglist);
+      return fold_builtin_tan (arglist, type);

     CASE_FLT_FN (BUILT_IN_ATAN):
       return fold_builtin_atan (arglist, type);
@@ -11170,3 +11171,53 @@ init_target_chars (void)
     }
   return true;
 }
+
+/* If argument ARG is a REAL_CST, call the one-argument mpfr function
+   FUNC on it and return the resulting value as a tree with type TYPE.
+   The mpfr precision is set to the precision of TYPE.  We assume that
+   function FUNC returns zero if the result could be calculated
+   exactly within the requested precision.  */
+
+static tree
+do_mpfr_arg1 (tree arg, tree type, int (*func)(mpfr_ptr, mpfr_srcptr, mp_rnd_t))
+{
+  tree result = NULL_TREE;
+
+  STRIP_NOPS (arg);
+
+  if (TREE_CODE (arg) == REAL_CST && ! TREE_CONSTANT_OVERFLOW (arg))
+    {
+      REAL_VALUE_TYPE r = TREE_REAL_CST (arg);
+
+      if (!real_isnan (&r) && !real_isinf (&r))
+        {
+	  const enum machine_mode mode = TYPE_MODE (type);
+	  const int prec = REAL_MODE_FORMAT (mode)->p;
+	  int exact;
+	  mpfr_t m;
+
+	  /* Initialize to double the requested precision.  */
+	  mpfr_init2 (m, prec*2);
+
+	  mpfr_from_real (m, &r);
+	  exact = func (m, m, GMP_RNDN);
+
+	  /* Now round to the requested precision.  */
+	  exact |= mpfr_prec_round (m, prec, GMP_RNDN);
+
+	  /* Proceed iff we get a normal number, i.e. not NaN or Inf.
+	     If -frounding-math is set, proceed iff the result of
+	     calling FUNC was exact, i.e. FUNC returned zero.  */
+	  if (mpfr_number_p (m)
+	      && (! flag_rounding_math || exact == 0))
+	  {
+	    real_from_mpfr (&r, m);
+	    real_convert (&r, mode, &r);
+	    result = build_real (type, r);
+	  }
+	  mpfr_clear (m);
+	}
+    }
+
+  return result;
+}
diff -rup orig/egcc-SVN20061016/gcc/real.c egcc-SVN20061016/gcc/real.c
--- orig/egcc-SVN20061016/gcc/real.c	2006-03-16 20:01:37.000000000 -0500
+++ egcc-SVN20061016/gcc/real.c	2006-10-19 20:44:30.534629192 -0400
@@ -4922,3 +4922,47 @@ real_copysign (REAL_VALUE_TYPE *r, const
   r->sign = x->sign;
 }

+/* Convert from REAL_VALUE_TYPE to MPFR.  The caller is responsible
+   for initializing and clearing the MPFR parmeter.  */
+
+void
+mpfr_from_real (mpfr_ptr m, const REAL_VALUE_TYPE *r)
+{
+  /* We use a string as an intermediate type.  */
+  char buf[128];
+
+  real_to_hexadecimal (buf, r, sizeof (buf), 0, 1);
+  /* mpfr_set_str() parses hexadecimal floats from strings in the same
+     format that GCC will output them.  Nothing extra is needed.  */
+  gcc_assert (mpfr_set_str (m, buf, 16, GMP_RNDN) == 0);
+}
+
+/* Convert from MPFR to REAL_VALUE_TYPE.  */
+
+void
+real_from_mpfr (REAL_VALUE_TYPE *r, mpfr_srcptr m)
+{
+  /* We use a string as an intermediate type.  */
+  char buf[128], *rstr;
+  mp_exp_t exp;
+
+  rstr = mpfr_get_str (NULL, &exp, 16, 0, m, GMP_RNDN);
+
+  /* The additional 12 chars add space for the sprintf below.  This
+     leaves 6 digits for the exponent which is supposedly enough.  */
+  gcc_assert (rstr != NULL && strlen (rstr) < sizeof (buf) - 12);
+
+  /* REAL_VALUE_ATOF expects the exponent for mantissa * 2**exp,
+     mpfr_get_str returns the exponent for mantissa * 16**exp, adjust
+     for that.  */
+  exp *= 4;
+
+  if (rstr[0] == '-')
+    sprintf (buf, "-0x.%sp%d", &rstr[1], (int) exp);
+  else
+    sprintf (buf, "0x.%sp%d", rstr, (int) exp);
+
+  mpfr_free_str (rstr);
+
+  real_from_string (r, buf);
+}
diff -rup orig/egcc-SVN20061016/gcc/real.h egcc-SVN20061016/gcc/real.h
--- orig/egcc-SVN20061016/gcc/real.h	2006-01-23 00:24:02.000000000 -0500
+++ egcc-SVN20061016/gcc/real.h	2006-10-19 13:06:22.016639937 -0400
@@ -22,6 +22,8 @@
 #ifndef GCC_REAL_H
 #define GCC_REAL_H

+#include <gmp.h>
+#include <mpfr.h>
 #include "machmode.h"

 /* An expanded form of the represented number.  */
@@ -425,4 +427,10 @@ extern void real_round (REAL_VALUE_TYPE
 /* Set the sign of R to the sign of X.  */
 extern void real_copysign (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);

+/* Convert between MPFR and REAL_VALUE_TYPE.  The caller is
+   responsible for initializing and clearing the MPFR parameter.  */
+
+extern void real_from_mpfr (REAL_VALUE_TYPE *, mpfr_srcptr);
+extern void mpfr_from_real (mpfr_ptr, const REAL_VALUE_TYPE *);
+
 #endif /* ! GCC_REAL_H */