re PR fortran/36117 (Use MPFR for bessel function (optimization, rejects valid F2008))

2008-05-06  Tobias Burnus  <burnus@net-b.de>

        PR fortran/36117
        * intrinsic.c (add_functions): Call gfc_simplify_bessel_*.
        * intrinsic.h: Add prototypes for gfc_simplify_bessel_*.
        * simplify.c (gfc_simplify_bessel_j0,gfc_simplify_bessel_j1,
        gfc_simplify_bessel_jn,gfc_simplify_bessel_y0,
        gfc_simplify_bessel_y1,gfc_simplify_bessel_yn): New.

2008-05-06  Tobias Burnus  <burnus@net-b.de>

        PR fortran/36117
        * gfortran.dg/bessel_2.f90: New.

From-SVN: r134988

gcc/fortran/ChangeLog

+2008-05-06  Tobias Burnus  <burnus@net-b.de>
+
+	PR fortran/36117
+	* intrinsic.c (add_functions): Call gfc_simplify_bessel_*.
+	* intrinsic.h: Add prototypes for gfc_simplify_bessel_*.
+	* simplify.c (gfc_simplify_bessel_j0,gfc_simplify_bessel_j1,
+	gfc_simplify_bessel_jn,gfc_simplify_bessel_y0,
+	gfc_simplify_bessel_y1,gfc_simplify_bessel_yn): New.
+
 2008-05-03  Janus Weil  <jaydub66@gmail.com>
 
 	* misc.c (gfc_clear_ts): Set interface to NULL.

gcc/fortran/intrinsic.c

@@ -1095,73 +1095,73 @@ add_functions (void)
   
   /* Bessel and Neumann functions for G77 compatibility.  */
   add_sym_1 ("besj0", GFC_ISYM_J0, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
-	     gfc_check_fn_r, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_r, gfc_simplify_bessel_j0, gfc_resolve_g77_math1,
 	     x, BT_REAL, dr, REQUIRED);
 
   make_alias ("bessel_j0", GFC_STD_F2008);
 
   add_sym_1 ("dbesj0", GFC_ISYM_J0, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dd, GFC_STD_GNU,
-	     gfc_check_fn_d, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_d, gfc_simplify_bessel_j0, gfc_resolve_g77_math1,
 	     x, BT_REAL, dd, REQUIRED);
 
   make_generic ("bessel_j0", GFC_ISYM_J0, GFC_STD_F2008);
 
   add_sym_1 ("besj1", GFC_ISYM_J1, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
-	     gfc_check_fn_r, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_r, gfc_simplify_bessel_j1, gfc_resolve_g77_math1,
 	     x, BT_REAL, dr, REQUIRED);
 
   make_alias ("bessel_j1", GFC_STD_F2008);
 
   add_sym_1 ("dbesj1", GFC_ISYM_J1, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dd, GFC_STD_GNU,
-	     gfc_check_fn_d, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_d, gfc_simplify_bessel_j1, gfc_resolve_g77_math1,
 	     x, BT_REAL, dd, REQUIRED);
 
   make_generic ("bessel_j1", GFC_ISYM_J1, GFC_STD_F2008);
 
   add_sym_2 ("besjn", GFC_ISYM_JN, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
-	     gfc_check_besn, NULL, gfc_resolve_besn,
+	     gfc_check_besn, gfc_simplify_bessel_jn, gfc_resolve_besn,
 	     n, BT_INTEGER, di, REQUIRED, x, BT_REAL, dr, REQUIRED);
 
   make_alias ("bessel_jn", GFC_STD_F2008);
 
   add_sym_2 ("dbesjn", GFC_ISYM_JN, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dd, GFC_STD_GNU,
-	     gfc_check_besn, NULL, gfc_resolve_besn,
+	     gfc_check_besn, gfc_simplify_bessel_jn, gfc_resolve_besn,
 	     n, BT_INTEGER, di, REQUIRED, x, BT_REAL, dd, REQUIRED);
 
   make_generic ("bessel_jn", GFC_ISYM_JN, GFC_STD_F2008);
 
   add_sym_1 ("besy0", GFC_ISYM_Y0, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
-	     gfc_check_fn_r, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_r, gfc_simplify_bessel_y0, gfc_resolve_g77_math1,
 	     x, BT_REAL, dr, REQUIRED);
 
   make_alias ("bessel_y0", GFC_STD_F2008);
 
   add_sym_1 ("dbesy0", GFC_ISYM_Y0, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dd, GFC_STD_GNU,
-	     gfc_check_fn_d, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_d, gfc_simplify_bessel_y0, gfc_resolve_g77_math1,
 	     x, BT_REAL, dd, REQUIRED);
 
   make_generic ("bessel_y0", GFC_ISYM_Y0, GFC_STD_F2008);
 
   add_sym_1 ("besy1", GFC_ISYM_Y1, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
-	     gfc_check_fn_r, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_r, gfc_simplify_bessel_y1, gfc_resolve_g77_math1,
 	     x, BT_REAL, dr, REQUIRED);
 
   make_alias ("bessel_y1", GFC_STD_F2008);
 
   add_sym_1 ("dbesy1", GFC_ISYM_Y1, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dd, GFC_STD_GNU,
-	     gfc_check_fn_d, NULL, gfc_resolve_g77_math1,
+	     gfc_check_fn_d, gfc_simplify_bessel_y1, gfc_resolve_g77_math1,
 	     x, BT_REAL, dd, REQUIRED);
 
   make_generic ("bessel_y1", GFC_ISYM_Y1, GFC_STD_F2008);
 
   add_sym_2 ("besyn", GFC_ISYM_YN, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
-	     gfc_check_besn, NULL, gfc_resolve_besn,
+	     gfc_check_besn, gfc_simplify_bessel_yn, gfc_resolve_besn,
 	     n, BT_INTEGER, di, REQUIRED, x, BT_REAL, dr, REQUIRED);
 
   make_alias ("bessel_yn", GFC_STD_F2008);
 
   add_sym_2 ("dbesyn", GFC_ISYM_YN, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dd, GFC_STD_GNU,
-	     gfc_check_besn, NULL, gfc_resolve_besn,
+	     gfc_check_besn, gfc_simplify_bessel_yn, gfc_resolve_besn,
 	     n, BT_INTEGER, di, REQUIRED, x, BT_REAL, dd, REQUIRED);
 
   make_generic ("bessel_yn", GFC_ISYM_YN, GFC_STD_F2008);

gcc/fortran/intrinsic.h

@@ -208,6 +208,12 @@ gfc_expr *gfc_simplify_asinh (gfc_expr *);
 gfc_expr *gfc_simplify_atan (gfc_expr *);
 gfc_expr *gfc_simplify_atanh (gfc_expr *);
 gfc_expr *gfc_simplify_atan2 (gfc_expr *, gfc_expr *);
+gfc_expr *gfc_simplify_bessel_j0 (gfc_expr *);
+gfc_expr *gfc_simplify_bessel_j1 (gfc_expr *);
+gfc_expr *gfc_simplify_bessel_jn (gfc_expr *, gfc_expr *);
+gfc_expr *gfc_simplify_bessel_y0 (gfc_expr *);
+gfc_expr *gfc_simplify_bessel_y1 (gfc_expr *);
+gfc_expr *gfc_simplify_bessel_yn (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_bit_size (gfc_expr *);
 gfc_expr *gfc_simplify_btest (gfc_expr *, gfc_expr *);
 gfc_expr *gfc_simplify_ceiling (gfc_expr *, gfc_expr *);

gcc/fortran/simplify.c

@@ -637,6 +637,130 @@ gfc_simplify_atan2 (gfc_expr *y, gfc_expr *x)
 
 
 gfc_expr *
+gfc_simplify_bessel_j0 (gfc_expr *x)
+{
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_j0 (result->value.real, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "BESSEL_J0");
+#else
+  return NULL;
+#endif
+}
+
+
+gfc_expr *
+gfc_simplify_bessel_j1 (gfc_expr *x)
+{
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_j1 (result->value.real, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "BESSEL_J1");
+#else
+  return NULL;
+#endif
+}
+
+
+gfc_expr *
+gfc_simplify_bessel_jn (gfc_expr *order, gfc_expr *x)
+{
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+  gfc_expr *result;
+  long n;
+
+  if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  n = mpz_get_si (order->value.integer);
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_jn (result->value.real, n, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "BESSEL_JN");
+#else
+  return NULL;
+#endif
+}
+
+
+gfc_expr *
+gfc_simplify_bessel_y0 (gfc_expr *x)
+{
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_y0 (result->value.real, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "BESSEL_Y0");
+#else
+  return NULL;
+#endif
+}
+
+
+gfc_expr *
+gfc_simplify_bessel_y1 (gfc_expr *x)
+{
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+  gfc_expr *result;
+
+  if (x->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_y1 (result->value.real, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "BESSEL_Y1");
+#else
+  return NULL;
+#endif
+}
+
+
+gfc_expr *
+gfc_simplify_bessel_yn (gfc_expr *order, gfc_expr *x)
+{
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+  gfc_expr *result;
+  long n;
+
+  if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
+    return NULL;
+
+  n = mpz_get_si (order->value.integer);
+  result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+  gfc_set_model_kind (x->ts.kind);
+  mpfr_yn (result->value.real, n, x->value.real, GFC_RND_MODE);
+
+  return range_check (result, "BESSEL_YN");
+#else
+  return NULL;
+#endif
+}
+
+
+gfc_expr *
 gfc_simplify_bit_size (gfc_expr *e)
 {
   gfc_expr *result;

gcc/testsuite/ChangeLog

+2008-05-06  Tobias Burnus  <burnus@net-b.de>
+
+	PR fortran/36117
+	* gfortran.dg/bessel_2.f90: New.
+
 2008-05-06  Olivier Hainque  <hainque@adacore.com>
 
 	* gnat.dg/fatp_sra.adb: New test.

gcc/testsuite/gfortran.dg/bessel_2.f90

+! { dg-do compile }
+! PR fortran/36117
+!
+! This program will fail for MPFR < 2.3.0
+!
+! Based on a test by James Van Buskirk.
+!
+program bug3
+   implicit none
+   real, parameter :: Qarg1 = 1.7
+   integer, parameter :: k2 = kind(BESJ0(Qarg1))
+   integer, parameter :: is_int = 1-1/(2+0*BESJ0(Qarg1))*2
+   integer, parameter :: kind_if_real = &
+      (1-is_int)*k2+is_int*kind(1.0)
+   complex :: z = cmplx(0,1,kind_if_real) ! FAILS
+   if (kind_if_real /= kind(Qarg1)) call abort ()
+end program bug3