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]: PR middle-end/30250, use MPFR for lgamma_r/gamma_r


This patch partly addresses PR middle-end/30250.  It uses MPFR to
evaluate the reentrant lgamma_r/gamma_r builtins at compile-time when
supplied with constant args.  The reentrant versions supply signgam as
a pointer parameter, and they were slightly easier to implement than
doing the regular lgamma functions, so I did them first.  (The plain
lgamma/gamma patch will come hopefully soon.)

The patch requires several others to apply cleanly:
http://gcc.gnu.org/ml/gcc-patches/2007-04/msg01624.html
http://gcc.gnu.org/ml/gcc-patches/2007-04/msg01663.html
http://gcc.gnu.org/ml/gcc-patches/2007-05/msg00297.html
http://gcc.gnu.org/ml/gcc-patches/2007-05/msg00320.html

Patch below tested on sparc-sun-solaris2.10 in conjuction with the
above four patches, no regressions and the new testcases all pass.  I
also ran a preliminary version of my plain lgamma/gamma patch against
the FP accuracy testsuite to check whether MPFR's lgamma function was
producing correct results back to GCC.  It was, all inputs generated
"perfect" results.

Okay for mainline?

		Thanks,
		--Kaveh

:ADDPATCH middle-end:


2007-05-07  Kaveh R. Ghazi  <ghazi@caip.rutgers.edu>

	PR middle-end/30250
	* builtins.c (do_mpfr_lgamma_r): New.
	(fold_builtin_2): Handle builtin gamma_r/lgamma_r.
	* tree.h (CASE_FLT_FN_REENT): New.

testsuite:
	* gcc.dg/torture/builtin-math-2.c: Add gamma_r/lgamma_r tests.
	* gcc.dg/torture/builtin-math-4.c: Likewise.

diff -rup orig/egcc-SVN20070505/gcc/builtins.c egcc-SVN20070505/gcc/builtins.c
--- orig/egcc-SVN20070505/gcc/builtins.c	2007-05-06 23:43:01.621873788 -0400
+++ egcc-SVN20070505/gcc/builtins.c	2007-05-07 00:55:32.966335476 -0400
@@ -234,6 +234,7 @@ static tree do_mpfr_bessel_n (tree, tree
 			      int (*)(mpfr_ptr, long, mpfr_srcptr, mp_rnd_t),
 			      const REAL_VALUE_TYPE *, bool);
 static tree do_mpfr_remquo (tree, tree, tree);
+static tree do_mpfr_lgamma_r (tree, tree, tree);
 #endif

 /* Return true if NODE should be considered for inline expansion regardless
@@ -9862,6 +9863,13 @@ fold_builtin_2 (tree fndecl, tree arg0,
           && validate_arg(arg1, REAL_TYPE))
         return do_mpfr_arg2 (arg0, arg1, type, mpfr_remainder);
     break;
+
+    CASE_FLT_FN_REENT (BUILT_IN_GAMMA): /* GAMMA_R */
+    CASE_FLT_FN_REENT (BUILT_IN_LGAMMA): /* LGAMMA_R */
+      if (validate_arg (arg0, REAL_TYPE)
+	  && validate_arg(arg1, POINTER_TYPE))
+	return do_mpfr_lgamma_r (arg0, arg1, type);
+    break;
 #endif

     CASE_FLT_FN (BUILT_IN_ATAN2):
@@ -12616,4 +12624,66 @@ do_mpfr_remquo (tree arg0, tree arg1, tr
     }
   return result;
 }
+
+/* If ARG is a REAL_CST, call mpfr_lgamma() 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 this mpfr function
+   returns zero if the result could be calculated exactly within the
+   requested precision.  In addition, the integer pointer represented
+   by ARG_SG will be dereferenced and set to the appropriate signgam
+   (-1,1) value.  */
+
+static tree
+do_mpfr_lgamma_r (tree arg, tree arg_sg, tree type)
+{
+  tree result = NULL_TREE;
+
+  STRIP_NOPS (arg);
+
+  /* To proceed, MPFR must exactly represent the target floating point
+     format, which only happens when the target base equals two.  */
+  if (REAL_MODE_FORMAT (TYPE_MODE (type))->b == 2
+      && TREE_CODE (arg) == REAL_CST && !TREE_OVERFLOW (arg))
+    {
+      const REAL_VALUE_TYPE *const ra = TREE_REAL_CST_PTR (arg);
+
+      /* In addition to NaN and Inf, the argument cannot be zero or a
+	 negative integer.  */
+      if (!real_isnan (ra) && !real_isinf (ra)
+	  && ra->cl != rvc_zero
+	  && !(real_isneg(ra) && real_isinteger(ra, TYPE_MODE (type))))
+        {
+	  const int prec = REAL_MODE_FORMAT (TYPE_MODE (type))->p;
+	  int inexact, sg;
+	  mpfr_t m;
+	  tree result_lg;
+
+	  mpfr_init2 (m, prec);
+	  mpfr_from_real (m, ra, GMP_RNDN);
+	  mpfr_clear_flags ();
+	  inexact = mpfr_lgamma (m, &sg, m, GMP_RNDN);
+	  result_lg = do_mpfr_ckconv (m, type, inexact);
+	  mpfr_clear (m);
+	  if (result_lg)
+	    {
+	      /* Dereference the arg_sg pointer argument.  */
+	      arg_sg = build_fold_indirect_ref (arg_sg);
+	      /* Proceed iff a valid pointer type was passed in.  */
+	      if (TYPE_MAIN_VARIANT (TREE_TYPE (arg_sg)) == integer_type_node)
+	        {
+		  /* Set the value. */
+		  tree result_sg = fold_build2 (MODIFY_EXPR,
+						TREE_TYPE (arg_sg), arg_sg,
+						 build_int_cst (NULL, sg));
+		  TREE_SIDE_EFFECTS (result_sg) = 1;
+		  /* Combine the signgam assignment with the lgamma result.  */
+		  result = non_lvalue (fold_build2 (COMPOUND_EXPR, type,
+						    result_sg, result_lg));
+		}
+	    }
+	}
+    }
+
+  return result;
+}
 #endif
diff -rup orig/egcc-SVN20070505/gcc/tree.h egcc-SVN20070505/gcc/tree.h
--- orig/egcc-SVN20070505/gcc/tree.h	2007-05-03 03:14:35.000000000 -0400
+++ egcc-SVN20070505/gcc/tree.h	2007-05-06 23:46:56.209230787 -0400
@@ -273,6 +273,7 @@ extern const char * built_in_names[(int)
 #define BUILTIN_ROOT_P(FN) (BUILTIN_SQRT_P (FN) || BUILTIN_CBRT_P (FN))

 #define CASE_FLT_FN(FN) case FN: case FN##F: case FN##L
+#define CASE_FLT_FN_REENT(FN) case FN##_R: case FN##F_R: case FN##L_R
 #define CASE_INT_FN(FN) case FN: case FN##L: case FN##LL

 /* An array of _DECL trees for the above.  */
diff -rup orig/egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-2.c egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-2.c
--- orig/egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-2.c	2007-05-06 23:42:44.833833597 -0400
+++ egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-2.c	2007-05-06 23:46:56.202474683 -0400
@@ -49,6 +49,13 @@ extern void fool (long double);
   fool (__builtin_remquol (ARG1##L, ARG2##L, &quo)); \
 } while (0)

+#define TESTIT_REENT(FUNC,ARG1) do { \
+  int sg; \
+  foof (__builtin_##FUNC##f_r (ARG1##F, &sg)); \
+  foo (__builtin_##FUNC##_r (ARG1, &sg)); \
+  fool (__builtin_##FUNC##l_r (ARG1##L, &sg)); \
+} while (0)
+
 void bar()
 {
   /* An argument of NaN is not evaluated at compile-time.  */
@@ -266,6 +273,21 @@ void bar()
   TESTIT2 (remainder, 1.0, -0.0);
   TESTIT2 (drem, 1.0, 0.0);
   TESTIT2 (drem, 1.0, -0.0);
+
+  /* The argument to lgamma* cannot be zero or a negative integer.  */
+  TESTIT_REENT (lgamma, -4.0); /* lgamma_r */
+  TESTIT_REENT (lgamma, -3.0); /* lgamma_r */
+  TESTIT_REENT (lgamma, -2.0); /* lgamma_r */
+  TESTIT_REENT (lgamma, -1.0); /* lgamma_r */
+  TESTIT_REENT (lgamma, -0.0); /* lgamma_r */
+  TESTIT_REENT (lgamma, 0.0); /* lgamma_r */
+
+  TESTIT_REENT (gamma, -4.0); /* gamma_r */
+  TESTIT_REENT (gamma, -3.0); /* gamma_r */
+  TESTIT_REENT (gamma, -2.0); /* gamma_r */
+  TESTIT_REENT (gamma, -1.0); /* gamma_r */
+  TESTIT_REENT (gamma, -0.0); /* gamma_r */
+  TESTIT_REENT (gamma, 0.0); /* gamma_r */
 }

 /* { dg-final { scan-tree-dump-times "exp2 " 9 "original" } } */
@@ -340,4 +362,10 @@ void bar()
 /* { dg-final { scan-tree-dump-times "drem " 2 "original" } } */
 /* { dg-final { scan-tree-dump-times "dremf" 2 "original" } } */
 /* { dg-final { scan-tree-dump-times "dreml" 2 "original" } } */
+/* { dg-final { scan-tree-dump-times "lgamma_r " 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "lgammaf_r" 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "lgammal_r" 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "_gamma_r " 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "_gammaf_r" 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "_gammal_r" 6 "original" } } */
 /* { dg-final { cleanup-tree-dump "original" } } */
diff -rup orig/egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-4.c egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-4.c
--- orig/egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-4.c	2007-05-06 23:42:44.835144977 -0400
+++ egcc-SVN20070505/gcc/testsuite/gcc.dg/torture/builtin-math-4.c	2007-05-07 00:52:46.941083027 -0400
@@ -101,6 +101,44 @@ extern void link_error(int);
     link_error(__LINE__); \
   } while (0)

+/* Test that FUNC(ARG,&SG) == (RES) && SG == RES_SG.  */
+#define TESTIT_LGAMMA_REENT(FUNC,ARG,RES,RES_SG) do { \
+  int sg; \
+  sg = 123; \
+  if (__builtin_##FUNC##f_r(ARG##F,&sg) != RES##F \
+      || sg != RES_SG \
+      || CKSGN_F(__builtin_##FUNC##f_r(ARG##F,&sg),RES##F)) \
+    link_error(__LINE__); \
+  sg = 123; \
+  if (__builtin_##FUNC##_r(ARG,&sg) != RES \
+      || sg != RES_SG \
+      || CKSGN(__builtin_##FUNC##_r(ARG,&sg),RES)) \
+    link_error(__LINE__); \
+  sg = 123; \
+  if (__builtin_##FUNC##l_r(ARG##L,&sg) != RES##L \
+      || sg != RES_SG \
+      || CKSGN_L(__builtin_##FUNC##l_r(ARG##L,&sg),RES##L)) \
+    link_error(__LINE__); \
+  } while (0)
+
+/* Range test, check that (LOW) < FUNC(ARG,&SG) < (HI), and also test
+   that SG == RES_SG.  */
+#define TESTIT_LGAMMA_REENT_R(FUNC,ARG,LOW,HI,RES_SG) do { \
+  int sg; \
+  sg = 123; \
+  if (__builtin_##FUNC##f_r(ARG,&sg) <= (LOW) || __builtin_##FUNC##f_r(ARG,&sg) >= (HI) \
+      || sg != RES_SG) \
+    link_error(__LINE__); \
+  sg = 123; \
+  if (__builtin_##FUNC##_r(ARG,&sg) <= (LOW) || __builtin_##FUNC##_r(ARG,&sg) >= (HI) \
+      || sg != RES_SG) \
+    link_error(__LINE__); \
+  sg = 123; \
+  if (__builtin_##FUNC##l_r(ARG,&sg) <= (LOW) || __builtin_##FUNC##l_r(ARG,&sg) >= (HI) \
+      || sg != RES_SG) \
+    link_error(__LINE__); \
+  } while (0)
+
 int main (void)
 {
 #ifdef __OPTIMIZE__
@@ -247,7 +285,27 @@ int main (void)
     MAXIT(remquol, -(__INT_MAX__+1.0L), 0);
     MAXIT(remquol, -(__INT_MAX__+2.0L), -1);
   }
-#endif

+  /* These tests rely on propagating the variable sg which contains
+     signgam.  This happens only when optimization is turned on.  */
+  TESTIT_LGAMMA_REENT_R (lgamma, -2.5, -0.06, -0.05, -1); /* lgamma_r(-2.5) == -0.056... */
+  TESTIT_LGAMMA_REENT_R (lgamma, -1.5, 0.86, 0.87, 1); /* lgamma_r(-1.5) == 0.860... */
+  TESTIT_LGAMMA_REENT_R (lgamma, -0.5, 1.26, 1.27, -1); /* lgamma_r(-0.5) == 1.265... */
+  TESTIT_LGAMMA_REENT_R (lgamma, 0.5, 0.57, 0.58, 1); /* lgamma_r(0.5) == 0.572... */
+  TESTIT_LGAMMA_REENT (lgamma, 1.0, 0.0, 1); /* lgamma_r(1) == 0 */
+  TESTIT_LGAMMA_REENT_R (lgamma, 1.5, -0.13, -0.12, 1); /* lgamma_r(1.5) == -0.120... */
+  TESTIT_LGAMMA_REENT (lgamma, 2.0, 0.0, 1); /* lgamma_r(2) == 0 */
+  TESTIT_LGAMMA_REENT_R (lgamma, 2.5, 0.28, 0.29, 1); /* lgamma_r(2.5) == 0.284... */
+
+  TESTIT_LGAMMA_REENT_R (gamma, -2.5, -0.06, -0.05, -1); /* gamma_r(-2.5) == -0.056... */
+  TESTIT_LGAMMA_REENT_R (gamma, -1.5, 0.86, 0.87, 1); /* gamma_r(-1.5) == 0.860... */
+  TESTIT_LGAMMA_REENT_R (gamma, -0.5, 1.26, 1.27, -1); /* gamma_r(-0.5) == 1.265... */
+  TESTIT_LGAMMA_REENT_R (gamma, 0.5, 0.57, 0.58, 1); /* gamma_r(0.5) == 0.572... */
+  TESTIT_LGAMMA_REENT (gamma, 1.0, 0.0, 1); /* gamma_r(1) == 0 */
+  TESTIT_LGAMMA_REENT_R (gamma, 1.5, -0.13, -0.12, 1); /* gamma_r(1.5) == -0.120... */
+  TESTIT_LGAMMA_REENT (gamma, 2.0, 0.0, 1); /* gamma_r(2) == 0 */
+  TESTIT_LGAMMA_REENT_R (gamma, 2.5, 0.28, 0.29, 1); /* gamma_r(2.5) == 0.284... */
+#endif
+
   return 0;
 }


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