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]: middle-end: Use MPFR for Bessel functions y0, y1, yn


This patch uses MPFR to evaluate at compile-time the y* Bessel functions.
This patch relies on the previous j* Bessel patch posted here:

http://gcc.gnu.org/ml/gcc-patches/2007-04/msg01624.html

Tested on sparc-sun-solaris2.10 with the mpfr-2.3.0 trunk, no regressions.

I also tested this against the famous FP accuracy suite I've been using
and it reported there were supposedly some problems with the results.
When I discussed this with the author of the MPFR Bessel routine, he
suspected it was a bug in Maple because Mathematica and GP/Pari both
agreed with MPFR's numbers.  (The FP accuracy testsuite creates the
"expected" values using Maple.)  After reporting this problem, I was told
the folks at Maplesoft acknowledged it was a bug in Maple's BesselY code
and resolved to fix it in a future Maple release.

So, okay for mainline?

		Thanks,
		--Kaveh


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

	* builtins.c (fold_builtin_1): Handle y0, y1.
	(fold_builtin_2): Handle yn.

testsuite:
	* gcc.dg/torture/builtin-math-2.c: Test y0, y1, yn.
	* gcc.dg/torture/builtin-math-4.c: Likewise.

diff -rup orig/egcc-SVN20070422/gcc/builtins.c egcc-SVN20070422/gcc/builtins.c
--- orig/egcc-SVN20070422/gcc/builtins.c	2007-04-24 15:57:02.177741358 -0400
+++ egcc-SVN20070422/gcc/builtins.c	2007-04-24 16:47:30.046334146 -0400
@@ -9710,6 +9710,18 @@ fold_builtin_1 (tree fndecl, tree arg0,
 	return do_mpfr_arg1 (arg0, type, mpfr_j1,
 			     NULL, NULL, 0);
     break;
+
+    CASE_FLT_FN (BUILT_IN_Y0):
+      if (validate_arg (arg0, REAL_TYPE))
+	return do_mpfr_arg1 (arg0, type, mpfr_y0,
+			     &dconst0, NULL, false);
+    break;
+
+    CASE_FLT_FN (BUILT_IN_Y1):
+      if (validate_arg (arg0, REAL_TYPE))
+	return do_mpfr_arg1 (arg0, type, mpfr_y1,
+			     &dconst0, NULL, false);
+    break;
 #endif

     CASE_FLT_FN (BUILT_IN_NAN):
@@ -9828,6 +9840,13 @@ fold_builtin_2 (tree fndecl, tree arg0,
 	  && validate_arg (arg1, REAL_TYPE))
 	return do_mpfr_bessel_n (arg0, arg1, type, mpfr_jn, NULL, 0);
     break;
+
+    CASE_FLT_FN (BUILT_IN_YN):
+      if (validate_arg (arg0, INTEGER_TYPE)
+	  && validate_arg (arg1, REAL_TYPE))
+	return do_mpfr_bessel_n (arg0, arg1, type, mpfr_yn,
+				 &dconst0, false);
+    break;
 #endif

     CASE_FLT_FN (BUILT_IN_ATAN2):
diff -rup orig/egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-2.c egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-2.c
--- orig/egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-2.c	2007-02-23 09:49:41.000000000 -0500
+++ egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-2.c	2007-04-24 16:30:25.320364610 -0400
@@ -24,6 +24,12 @@ extern void fool (long double);
   fool (__builtin_##FUNC##l (ARG1##L, ARG2##L)); \
 } while (0)

+#define TESTIT2_I1(FUNC, ARG1, ARG2) do { \
+  foof (__builtin_##FUNC##f (ARG1, ARG2##F)); \
+  foo (__builtin_##FUNC (ARG1, ARG2)); \
+  fool (__builtin_##FUNC##l (ARG1, ARG2##L)); \
+} while (0)
+
 #define TESTIT2_I2ALL(FUNC, ARGF, MAXF, ARGD, MAXD, ARGLD, MAXLD) do { \
   foof (__builtin_##FUNC##f (ARGF, MAXF)); \
   foo (__builtin_##FUNC (ARGD, MAXD)); \
@@ -228,6 +234,24 @@ void bar()
   foof (__builtin_ilogbf (-__builtin_nanf("")));
   foo (__builtin_ilogb (-__builtin_nan("")));
   fool (__builtin_ilogbl (-__builtin_nanl("")));
+
+  /* The y* arg must be [0 ... Inf] EXclusive.  */
+  TESTIT (y0, -1.0);
+  TESTIT (y0, 0.0);
+  TESTIT (y0, -0.0);
+
+  TESTIT (y1, -1.0);
+  TESTIT (y1, 0.0);
+  TESTIT (y1, -0.0);
+
+  TESTIT2_I1 (yn, 2, -1.0);
+  TESTIT2_I1 (yn, 2, 0.0);
+  TESTIT2_I1 (yn, 2, -0.0);
+
+  TESTIT2_I1 (yn, -3, -1.0);
+  TESTIT2_I1 (yn, -3, 0.0);
+  TESTIT2_I1 (yn, -3, -0.0);
+
 }

 /* { dg-final { scan-tree-dump-times "exp2 " 9 "original" } } */
@@ -284,4 +308,13 @@ void bar()
 /* { dg-final { scan-tree-dump-times "ilogb " 6 "original" } } */
 /* { dg-final { scan-tree-dump-times "ilogbf" 6 "original" } } */
 /* { dg-final { scan-tree-dump-times "ilogbl" 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "y0 " 3 "original" } } */
+/* { dg-final { scan-tree-dump-times "y0f" 3 "original" } } */
+/* { dg-final { scan-tree-dump-times "y0l" 3 "original" } } */
+/* { dg-final { scan-tree-dump-times "y1 " 3 "original" } } */
+/* { dg-final { scan-tree-dump-times "y1f" 3 "original" } } */
+/* { dg-final { scan-tree-dump-times "y1l" 3 "original" } } */
+/* { dg-final { scan-tree-dump-times "yn " 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "ynf" 6 "original" } } */
+/* { dg-final { scan-tree-dump-times "ynl" 6 "original" } } */
 /* { dg-final { cleanup-tree-dump "original" } } */
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-24 15:57:02.179998752 -0400
+++ egcc-SVN20070422/gcc/testsuite/gcc.dg/torture/builtin-math-4.c	2007-04-24 17:03:46.824517068 -0400
@@ -111,5 +111,24 @@ int main (void)
   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... */

+  TESTIT_R (y0, 5.0, -0.31, -0.30); /* y0(5) == -0.308... */
+  TESTIT_R (y0, 0.1, -1.54, -1.53); /* y0(0.1) == -1.534... */
+
+  TESTIT_R (y1, 5.0, 0.14, 0.15); /* y1(5) == 0.147... */
+  TESTIT_R (y1, 0.1, -6.46, -6.45); /* y1(0.1) == -6.458... */
+
+  TESTIT2_R (yn, -1, 3.0, -0.33, -0.32); /* yn(-1,3) == -0.324... */
+  TESTIT2_R (yn, -1, 0.25, 2.70, 2.71); /* yn(-1,0.25) == 2.704... */
+
+  TESTIT2_R (yn, 2, 4.0, 0.21, 0.22); /* yn(2,4) == 0.215... */
+  TESTIT2_R (yn, 2, 0.9, -1.95, -1.94); /* yn(2,0.9) == -1.945... */
+  TESTIT2_R (yn, -2, 4.0, 0.21, 0.22); /* yn(-2,4) == 0.215... */
+  TESTIT2_R (yn, -2, 0.9, -1.95, -1.94); /* yn(-2,0.9) == -1.945... */
+
+  TESTIT2_R (yn, 3, 6.0, 0.32, 0.33); /* yn(3,6) == 0.328... */
+  TESTIT2_R (yn, 3, 0.89, -8.03, -8.02); /* yn(3,0.89) == -8.020... */
+  TESTIT2_R (yn, -3, 8.0, -0.03, -0.02); /* yn(-3,8) == -0.026... */
+  TESTIT2_R (yn, -3, 0.99, 5.98, 5.99); /* yn(-3,0.99) == 5.982... */
+
   return 0;
 }


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