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]: Use MPFR for Bessel functions j0, j1 and jn


The patch below uses MPFR to enable compile-time evaluation for Bessel
functions j0, j1 and jn.  This functionality appears in the upcoming
mpfr-2.3.0.  I've wrapped the new code in version checks so we don't
have to change GCC's configure tests for minimum required MPFR.  GCC
will continue to work fine (simply without this functionality) if used
with older MPFR versions.  I expect mpfr-2.3.0 to be released long
before gcc-4.3 is, and I'll track both until the new mpfr it released
to ensure nothing breaks.

Patch tested on sparc-sun-solaris2.10 using mpfr-2.2.1 and again
separately with MPFR's svn trunk.  There were no regressions in either
case and the new testcase passes when used with the MPFR trunk.

I also tested it against the FP accuracy testsuite which has tests for
j0/j1 and all of the cases got "perfect" results.

Okay for mainline?

		Thanks,
		--Kaveh


2007-04-23  Kaveh R. Ghazi  <ghazi@caip.rutgers.edu>

	* builtins.c (do_mpfr_bessel_n): New.
	(fold_builtin_1): Handle BUILT_IN_J0 and BUILT_IN_J1.
	(fold_builtin_2): Handle BUILT_IN_JN.

testsuite:
	* gcc.dg/torture/builtin-math-4.c: New test.

diff -rup orig/egcc-SVN20070422/gcc/builtins.c egcc-SVN20070422/gcc/builtins.c
--- orig/egcc-SVN20070422/gcc/builtins.c	2007-04-05 20:02:19.000000000 -0400
+++ egcc-SVN20070422/gcc/builtins.c	2007-04-23 01:05:26.406280618 -0400
@@ -229,6 +229,11 @@ static tree do_mpfr_arg2 (tree, tree, tr
 static tree do_mpfr_arg3 (tree, tree, tree, tree,
 			  int (*)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t));
 static tree do_mpfr_sincos (tree, tree, tree);
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+static tree do_mpfr_bessel_n (tree, tree, tree,
+			      int (*)(mpfr_ptr, long, mpfr_srcptr, mp_rnd_t),
+			      const REAL_VALUE_TYPE *, bool);
+#endif

 /* Return true if NODE should be considered for inline expansion regardless
    of the optimization level.  This means whenever a function is invoked with
@@ -9693,6 +9698,20 @@ fold_builtin_1 (tree fndecl, tree arg0,
 			     &dconstm1, NULL, false);
     break;

+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+    CASE_FLT_FN (BUILT_IN_J0):
+      if (validate_arg (arg0, REAL_TYPE))
+	return do_mpfr_arg1 (arg0, type, mpfr_j0,
+			     NULL, NULL, 0);
+    break;
+
+    CASE_FLT_FN (BUILT_IN_J1):
+      if (validate_arg (arg0, REAL_TYPE))
+	return do_mpfr_arg1 (arg0, type, mpfr_j1,
+			     NULL, NULL, 0);
+    break;
+#endif
+
     CASE_FLT_FN (BUILT_IN_NAN):
     case BUILT_IN_NAND32:
     case BUILT_IN_NAND64:
@@ -9803,6 +9822,13 @@ fold_builtin_2 (tree fndecl, tree arg0,

   switch (fcode)
     {
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+    CASE_FLT_FN (BUILT_IN_JN):
+      if (validate_arg (arg0, INTEGER_TYPE)
+	  && validate_arg (arg1, REAL_TYPE))
+	return do_mpfr_bessel_n (arg0, arg1, type, mpfr_jn, NULL, 0);
+    break;
+#endif

     CASE_FLT_FN (BUILT_IN_ATAN2):
       if (validate_arg (arg0, REAL_TYPE)
@@ -12429,3 +12455,50 @@ do_mpfr_sincos (tree arg, tree arg_sinp,
     }
   return result;
 }
+
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+/* If argument ARG1 is an INTEGER_CST and ARG2 is a REAL_CST, call the
+   two-argument mpfr order N Bessel function FUNC on them 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_bessel_n (tree arg1, tree arg2, tree type,
+		  int (*func)(mpfr_ptr, long, mpfr_srcptr, mp_rnd_t),
+		  const REAL_VALUE_TYPE *min, bool inclusive)
+{
+  tree result = NULL_TREE;
+
+  STRIP_NOPS (arg1);
+  STRIP_NOPS (arg2);
+
+  /* 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
+      && host_integerp (arg1, 0)
+      && TREE_CODE (arg2) == REAL_CST && !TREE_OVERFLOW (arg2))
+    {
+      const HOST_WIDE_INT n = tree_low_cst(arg1, 0);
+      const REAL_VALUE_TYPE *const ra = &TREE_REAL_CST (arg2);
+
+      if (n == (long)n
+	  && !real_isnan (ra) && !real_isinf (ra)
+	  && (!min || real_compare (inclusive ? GE_EXPR: GT_EXPR , ra, min)))
+        {
+	  const int prec = REAL_MODE_FORMAT (TYPE_MODE (type))->p;
+	  int inexact;
+	  mpfr_t m;
+
+	  mpfr_init2 (m, prec);
+	  mpfr_from_real (m, ra);
+	  mpfr_clear_flags ();
+	  inexact = func (m, n, m, GMP_RNDN);
+	  result = do_mpfr_ckconv (m, type, inexact);
+	  mpfr_clear (m);
+	}
+    }
+
+  return result;
+}
+#endif
diff -rup orig/egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-4.c egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-4.c
--- orig/egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-4.c	2007-04-23 01:06:10.767708035 -0400
+++ egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-4.c	2007-04-23 03:22:07.222334148 -0400
@@ -0,0 +1,115 @@
+/* Copyright (C) 2007  Free Software Foundation.
+
+   Verify that built-in math function constant folding of constant
+   arguments is correctly performed by the compiler.  This testcase is
+   for functionality that was available as of mpfr-2.3.0.
+
+   Origin: Kaveh R. Ghazi,  April 23, 2007.  */
+
+/* { dg-do link } */
+
+/* All references to link_error should go away at compile-time.  */
+extern void link_error(int);
+
+/* Return TRUE if the sign of X != sign of Y.  This is important when
+   comparing signed zeros.  */
+#define CKSGN_F(X,Y) \
+  (__builtin_copysignf(1.0F,(X)) != __builtin_copysignf(1.0F,(Y)))
+#define CKSGN(X,Y) \
+  (__builtin_copysign(1.0,(X)) != __builtin_copysign(1.0,(Y)))
+#define CKSGN_L(X,Y) \
+  (__builtin_copysignl(1.0L,(X)) != __builtin_copysignl(1.0L,(Y)))
+
+/* Test that FUNC(ARG) == (RES).  */
+#define TESTIT(FUNC,ARG,RES) do { \
+  if (__builtin_##FUNC##f(ARG##F) != RES##F \
+      || CKSGN_F(__builtin_##FUNC##f(ARG##F),RES##F)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC(ARG) != RES \
+      || CKSGN(__builtin_##FUNC(ARG),RES)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC##l(ARG##L) != RES##L \
+      || CKSGN_L(__builtin_##FUNC##l(ARG##L),RES##L)) \
+    link_error(__LINE__); \
+  } while (0)
+
+/* Range test, check that (LOW) < FUNC(ARG) < (HI).  */
+#define TESTIT_R(FUNC,ARG,LOW,HI) do { \
+  if (__builtin_##FUNC##f(ARG) <= (LOW) || __builtin_##FUNC##f(ARG) >= (HI)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC(ARG) <= (LOW) || __builtin_##FUNC(ARG) >= (HI)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC##l(ARG) <= (LOW) || __builtin_##FUNC##l(ARG) >= (HI)) \
+    link_error(__LINE__); \
+  } while (0)
+
+/* Test that FUNC(ARG1, ARG2) == (RES).  */
+#define TESTIT2(FUNC,ARG1,ARG2,RES) do { \
+  if (__builtin_##FUNC##f(ARG1, ARG2##F) != RES##F \
+      || CKSGN_F(__builtin_##FUNC##f(ARG1,ARG2##F),RES##F)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC(ARG1, ARG2) != RES \
+      || CKSGN(__builtin_##FUNC(ARG1,ARG2),RES)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC##l(ARG1, ARG2##L) != RES##L \
+      || CKSGN_L(__builtin_##FUNC##l(ARG1,ARG2##L),RES##L)) \
+    link_error(__LINE__); \
+  } while (0)
+
+/* Range test, check that (LOW) < FUNC(ARG1,ARG2) < (HI).  */
+#define TESTIT2_R(FUNC,ARG1,ARG2,LOW,HI) do { \
+  if (__builtin_##FUNC##f(ARG1, ARG2##F) <= (LOW) \
+      || __builtin_##FUNC##f(ARG1, ARG2##F) >= (HI)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC(ARG1, ARG2) <= (LOW) \
+      || __builtin_##FUNC(ARG1, ARG2) >= (HI)) \
+    link_error(__LINE__); \
+  if (__builtin_##FUNC##l(ARG1, ARG2##L) <= (LOW) \
+      || __builtin_##FUNC##l(ARG1, ARG2##L) >= (HI)) \
+    link_error(__LINE__); \
+  } while (0)
+
+int main (void)
+{
+  TESTIT (j0, 0.0, 1.0); /* j0(0) == 1 */
+  TESTIT (j0, -0.0, 1.0); /* j0(-0) == 1 */
+  TESTIT_R (j0, 1.0, 0.765, 0.766); /* j0(1) == 0.7651... */
+  TESTIT_R (j0, -1.0, 0.765, 0.766); /* j0(-1) == 0.7651... */
+
+  TESTIT (j1, 0.0, 0.0); /* j1(0) == 0 */
+  TESTIT (j1, -0.0, -0.0); /* j1(-0) == -0 */
+  TESTIT_R (j1, 1.0, 0.44, 0.45); /* j1(1) == 0.440... */
+  TESTIT_R (j1, -1.0, -0.45, -0.44); /* j1(-1) == -0.440... */
+
+  TESTIT2 (jn, 5, 0.0, 0.0); /* jn(5,0) == 0 */
+  TESTIT2 (jn, 5, -0.0, -0.0); /* jn(5,-0) == -0 */
+  TESTIT2 (jn, 6, 0.0, 0.0); /* jn(6,0) == 0 */
+  TESTIT2 (jn, 6, -0.0, 0.0); /* jn(6,-0) == 0 */
+
+  TESTIT2 (jn, -5, 0.0, -0.0); /* jn(-5,0) == -0 */
+  TESTIT2 (jn, -5, -0.0, 0.0); /* jn(-5,-0) == 0 */
+  TESTIT2 (jn, -6, 0.0, 0.0); /* jn(-6,0) == 0 */
+  TESTIT2 (jn, -6, -0.0, 0.0); /* jn(-6,-0) == 0 */
+
+  TESTIT2_R (jn, 2, 1.0, 0.11, 0.12); /* jn(2,1) == 0.114... */
+  TESTIT2_R (jn, 2, -1.0, 0.11, 0.12); /* jn(2,-1) == 0.114... */
+  TESTIT2_R (jn, 3, 5.0, 0.36, 0.37); /* jn(3,5) == 0.364... */
+  TESTIT2_R (jn, 3, -5.0, -0.37, -0.36); /* jn(3,-5) == -0.364... */
+
+  TESTIT2_R (jn, -2, 1.0, 0.11, 0.12); /* jn(-2,1) == 0.114... */
+  TESTIT2_R (jn, -2, -1.0, 0.11, 0.12); /* jn(-2,-1) == 0.114... */
+  TESTIT2_R (jn, -3, 5.0, -0.37, -0.36); /* jn(-3,5) == -0.364... */
+  TESTIT2_R (jn, -3, -5.0, 0.36, 0.37); /* jn(-3,-5) == 0.364... */
+
+  TESTIT2_R (jn, 4, 3.5, 0.20, 0.21); /* jn(4,3.5) == 0.204... */
+  TESTIT2_R (jn, 4, -3.5, 0.20, 0.21); /* jn(4,-3.5) == 0.204... */
+  TESTIT2_R (jn, 5, 4.6, 0.20, 0.21); /* jn(5,4.6) == 0.207... */
+  TESTIT2_R (jn, 5, -4.6, -0.21, -0.20); /* jn(5,-4.6) == -0.207... */
+
+  TESTIT2_R (jn, -4, 3.5, 0.20, 0.21); /* jn(-4,3.5) == 0.204... */
+  TESTIT2_R (jn, -4, -3.5, 0.20, 0.21); /* jn(-4,-3.5) == 0.204... */
+  TESTIT2_R (jn, -5, 4.6, -0.21, -0.20); /* jn(-5,4.6) == -0.207... */
+  TESTIT2_R (jn, -5, -4.6, 0.20, 0.21); /* jn(-5,-4.6) == 0.207... */
+
+  return 0;
+}


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