Patchwork [Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic

login
register
mail settings
Submitter Tobias Burnus
Date Aug. 17, 2010, 9:14 p.m.
Message ID <4C6AFBAF.2030009@net-b.de>
Download mbox | patch
Permalink /patch/61971/
State New
Headers show

Comments

Tobias Burnus - Aug. 17, 2010, 9:14 p.m.
Steve Kargl wrote:
> I'll look over the patch in more detail in a few hours.

Thanks! Here is an updated version. With the help of Daniel K, I have 
found a solution to the constructor issue; JN now works as it should. YN 
also works, though for some of the examples, the result differs quite a 
bit from the result of mpfr_yn. I do not know whether I have chosen 
particularly difficult values or whether the algorithm has problems. 
Please check. (Cf. bessel_5.f90 - the ULP values given are the 
differences which maximally occur [though 1 ULP could be 0 ULP in some 
cases].)

Build and currently regtesting on x86-64-linux.
OK for the trunk?

Tobias

Patch

2010-08-17  Tobias Burnus  <burnus@net-b.de>

	PR fortran/36158
	PR fortran/33197
	* check.c (gfc_check_bessel_n2): New function.
	* gfortran.h (gfc_isym_id): Add GFC_ISYM_JN2 and GFC_ISYM_YN2.
	* intrinsic.c (add_functions): Add transformational version
	of the Bessel_jn/yn intrinsics.
	* intrinsic.h (gfc_check_bessel_n2,gfc_simplify_bessel_jn2,
	gfc_simplify_bessel_yn2): New prototypes.
	* intrinsic.texi (Bessel_jn, Bessel_yn): Document
	transformational variant.
	* simplify.c (gfc_simplify_bessel_jn, gfc_simplify_bessel_yn):
	Check for negative order.
	(gfc_simplify_bessel_n2,gfc_simplify_bessel_jn2,
	gfc_simplify_bessel_yn2): New functions.

2010-08-17  Tobias Burnus  <burnus@net-b.de>

	PR fortran/36158
	PR fortran/33197
	* gfortran.dg/bessel_3.f90: New.
	* gfortran.dg/bessel_4.f90: New.
	* gfortran.dg/bessel_5.f90: New.

Index: gcc/fortran/intrinsic.c
===================================================================
--- gcc/fortran/intrinsic.c	(Revision 163318)
+++ gcc/fortran/intrinsic.c	(Arbeitskopie)
@@ -1317,6 +1317,11 @@  add_functions (void)
 	     gfc_check_besn, gfc_simplify_bessel_jn, gfc_resolve_besn,
 	     n, BT_INTEGER, di, REQUIRED, x, BT_REAL, dd, REQUIRED);
 
+  add_sym_3 ("bessel_jn", GFC_ISYM_JN2, CLASS_TRANSFORMATIONAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_F2008,
+	     gfc_check_bessel_n2, gfc_simplify_bessel_jn2, NULL,
+	     "n1", BT_INTEGER, di, REQUIRED,"n2", BT_INTEGER, di, REQUIRED,
+	     x, BT_REAL, dr, 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,
@@ -1353,6 +1358,11 @@  add_functions (void)
 	     gfc_check_besn, gfc_simplify_bessel_yn, gfc_resolve_besn,
 	     n, BT_INTEGER, di, REQUIRED, x, BT_REAL, dd, REQUIRED);
 
+  add_sym_3 ("bessel_yn", GFC_ISYM_YN2, CLASS_TRANSFORMATIONAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_F2008,
+	     gfc_check_bessel_n2, gfc_simplify_bessel_yn2, NULL,
+	     "n1", BT_INTEGER, di, REQUIRED,"n2", BT_INTEGER, di, REQUIRED,
+	      x, BT_REAL, dr, REQUIRED);
+
   make_generic ("bessel_yn", GFC_ISYM_YN, GFC_STD_F2008);
 
   add_sym_1 ("bit_size", GFC_ISYM_BIT_SIZE, CLASS_INQUIRY, ACTUAL_NO, BT_INTEGER, di, GFC_STD_F95,
Index: gcc/fortran/intrinsic.h
===================================================================
--- gcc/fortran/intrinsic.h	(Revision 163318)
+++ gcc/fortran/intrinsic.h	(Arbeitskopie)
@@ -40,6 +40,7 @@  gfc_try gfc_check_associated (gfc_expr *
 gfc_try gfc_check_atan_2 (gfc_expr *, gfc_expr *);
 gfc_try gfc_check_atan2 (gfc_expr *, gfc_expr *);
 gfc_try gfc_check_besn (gfc_expr *, gfc_expr *);
+gfc_try gfc_check_bessel_n2 (gfc_expr *, gfc_expr *, gfc_expr *);
 gfc_try gfc_check_bitfcn (gfc_expr *, gfc_expr *);
 gfc_try gfc_check_char (gfc_expr *, gfc_expr *);
 gfc_try gfc_check_chdir (gfc_expr *);
@@ -223,9 +224,11 @@  gfc_expr *gfc_simplify_atan2 (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_jn2 (gfc_expr *, 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_bessel_yn2 (gfc_expr *, 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 *);
Index: gcc/fortran/gfortran.h
===================================================================
--- gcc/fortran/gfortran.h	(Revision 163318)
+++ gcc/fortran/gfortran.h	(Arbeitskopie)
@@ -422,6 +422,7 @@  enum gfc_isym_id
   GFC_ISYM_J0,
   GFC_ISYM_J1,
   GFC_ISYM_JN,
+  GFC_ISYM_JN2,
   GFC_ISYM_KILL,
   GFC_ISYM_KIND,
   GFC_ISYM_LBOUND,
@@ -531,7 +532,8 @@  enum gfc_isym_id
   GFC_ISYM_XOR,
   GFC_ISYM_Y0,
   GFC_ISYM_Y1,
-  GFC_ISYM_YN
+  GFC_ISYM_YN,
+  GFC_ISYM_YN2
 };
 typedef enum gfc_isym_id gfc_isym_id;
 
Index: gcc/fortran/check.c
===================================================================
--- gcc/fortran/check.c	(Revision 163318)
+++ gcc/fortran/check.c	(Arbeitskopie)
@@ -884,11 +884,47 @@  gfc_check_besn (gfc_expr *n, gfc_expr *x
 {
   if (type_check (n, 0, BT_INTEGER) == FAILURE)
     return FAILURE;
+  if (n->expr_type == EXPR_CONSTANT)
+    {
+      int i;
+      gfc_extract_int (n, &i);
+      if (i < 0 && gfc_notify_std (GFC_STD_GNU, "Extension: Negative argument "
+				   "N at %L", &n->where) == FAILURE)
+	return FAILURE;
+    }
 
   if (type_check (x, 1, BT_REAL) == FAILURE)
     return FAILURE;
 
   return SUCCESS;
+}
+
+
+/* Transformational version of the Bessel JN and YN functions.  */
+
+gfc_try
+gfc_check_bessel_n2 (gfc_expr *n1, gfc_expr *n2, gfc_expr *x)
+{
+  if (type_check (n1, 0, BT_INTEGER) == FAILURE)
+    return FAILURE;
+  if (scalar_check (n1, 0) == FAILURE)
+    return FAILURE;
+  if (nonnegative_check("N1", n1) == FAILURE)
+    return FAILURE;
+
+  if (type_check (n2, 1, BT_INTEGER) == FAILURE)
+    return FAILURE;
+  if (scalar_check (n2, 1) == FAILURE)
+    return FAILURE;
+  if (nonnegative_check("N2", n2) == FAILURE)
+    return FAILURE;
+
+  if (type_check (x, 2, BT_REAL) == FAILURE)
+    return FAILURE;
+  if (scalar_check (x, 2) == FAILURE)
+    return FAILURE;
+
+  return SUCCESS;
 }
 
 
Index: gcc/fortran/intrinsic.texi
===================================================================
--- gcc/fortran/intrinsic.texi	(Revision 163318)
+++ gcc/fortran/intrinsic.texi	(Arbeitskopie)
@@ -1630,23 +1630,31 @@  end program test_besj1
 @item @emph{Description}:
 @code{BESSEL_JN(N, X)} computes the Bessel function of the first kind of
 order @var{N} of @var{X}. This function is available under the name
-@code{BESJN} as a GNU extension.
+@code{BESJN} as a GNU extension.  If @var{N} and @var{X} are arrays,
+their ranks and shapes shall conform.  
 
-If both arguments are arrays, their ranks and shapes shall conform.
+@code{BESSEL_JN(N1, N2, X)} returns an array with the Bessel functions
+of the first kind of the orders @var{N1} to @var{N2}.
 
 @item @emph{Standard}:
-Fortran 2008 and later
+Fortran 2008 and later, negative @var{N}, @var{N1}, @var{N2} are
+allowed as GNU extension
 
 @item @emph{Class}:
-Elemental function
+Elemental function, except for the tranformational function
+@code{BESSEL_JN(N1, N2, X)}
 
 @item @emph{Syntax}:
 @code{RESULT = BESSEL_JN(N, X)}
+@code{RESULT = BESSEL_JN(N1, N2, X)}
 
 @item @emph{Arguments}:
 @multitable @columnfractions .15 .70
 @item @var{N} @tab Shall be a scalar or an array of type  @code{INTEGER}.
+@item @var{N1} @tab Shall be a scalar of type  @code{INTEGER}.
+@item @var{N2} @tab Shall be a scalar of type  @code{INTEGER}.
 @item @var{X} @tab Shall be a scalar or an array of type  @code{REAL}.
+for @code{BESSEL_JN(N1, N2, X)} it shall be scalar.
 @end multitable
 
 @item @emph{Return value}:
@@ -1778,23 +1786,31 @@  end program test_besy1
 @item @emph{Description}:
 @code{BESSEL_YN(N, X)} computes the Bessel function of the second kind of
 order @var{N} of @var{X}. This function is available under the name
-@code{BESYN} as a GNU extension.
+@code{BESYN} as a GNU extension.  If @var{N} and @var{X} are arrays,
+their ranks and shapes shall conform.  
 
-If both arguments are arrays, their ranks and shapes shall conform.
+@code{BESSEL_YN(N1, N2, X)} returns an array with the Bessel functions
+of the first kind of the orders @var{N1} to @var{N2}.
 
 @item @emph{Standard}:
-Fortran 2008 and later
+Fortran 2008 and later, negative @var{N}, @var{N1}, @var{N2} are
+allowed as GNU extension
 
 @item @emph{Class}:
-Elemental function
+Elemental function, except for the tranformational function
+@code{BESSEL_YN(N1, N2, X)}
 
 @item @emph{Syntax}:
 @code{RESULT = BESSEL_YN(N, X)}
+@code{RESULT = BESSEL_YN(N1, N2, X)}
 
 @item @emph{Arguments}:
 @multitable @columnfractions .15 .70
-@item @var{N} @tab Shall be a scalar or an array of type  @code{INTEGER}.
-@item @var{X} @tab Shall be a scalar or an array of type  @code{REAL}.
+@item @var{N} @tab Shall be a scalar or an array of type  @code{INTEGER} .
+@item @var{N1} @tab Shall be a scalar of type  @code{INTEGER}.
+@item @var{N2} @tab Shall be a scalar of type  @code{INTEGER}.
+@item @var{X} @tab Shall be a scalar or an array of type  @code{REAL};
+for @code{BESSEL_JN(N1, N2, X)} it shall be scalar.
 @end multitable
 
 @item @emph{Return value}:
Index: gcc/fortran/simplify.c
===================================================================
--- gcc/fortran/simplify.c	(Revision 163318)
+++ gcc/fortran/simplify.c	(Arbeitskopie)
@@ -1196,6 +1196,183 @@  gfc_simplify_bessel_jn (gfc_expr *order,
 }
 
 
+/* Simplify transformational form of JN and YN.  */
+
+static gfc_expr *
+gfc_simplify_bessel_n2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x,
+			bool jn)
+{
+  gfc_expr *result;
+  gfc_expr *e;
+  long n1, n2;
+  int i;
+  mpfr_t x2rev, last1, last2;
+
+  if (x->expr_type != EXPR_CONSTANT || order1->expr_type != EXPR_CONSTANT
+      || order2->expr_type != EXPR_CONSTANT)
+    {
+      gfc_error ("Sorry, non-constant transformational Bessel function at %L"
+		   " not yet supported", &order2->where);
+      return &gfc_bad_expr;
+    }
+
+  n1 = mpz_get_si (order1->value.integer);
+  n2 = mpz_get_si (order2->value.integer);
+  result = gfc_get_array_expr (x->ts.type, x->ts.kind, &x->where);
+  result->rank = 1;
+  result->shape = gfc_get_shape (1);
+  mpz_init_set_ui (result->shape[0], MAX (n2-n1+1, 0));
+
+  if (n2 < n1)
+    return result;
+
+  /* Special case: x == 0; it is J0(0.0) == 1, JN(N > 0, 0.0) == 0; and
+     YN(N, 0.0) = -Inf.  */
+
+  if (mpfr_cmp_ui (x->value.real, 0.0) == 0)
+    {
+      if (!jn && gfc_option.flag_range_check)
+	{
+	  gfc_error ("Result of BESSEL_YN is -INF at %L", &result->where);
+ 	  gfc_free_expr (result);
+	  return &gfc_bad_expr;
+	}
+
+      if (jn && n1 == 0)
+	{
+	  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+	  mpfr_set_ui (e->value.real, 1.0, GFC_RND_MODE);
+	  gfc_constructor_append_expr (&result->value.constructor, e,
+				       &x->where);
+	  n1++;
+	}
+
+      for (i = n1; i <= n2; i++)
+	{
+	  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+	  if (jn)
+	    mpfr_set_ui (e->value.real, 0.0, GFC_RND_MODE);
+	  else
+            mpfr_set_inf (e->value.real, -1);
+	  gfc_constructor_append_expr (&result->value.constructor, e,
+				       &x->where);
+	}
+
+      return result;
+    }
+
+  /* Use the faster but more verbose recursion algorithm. Bessel functions
+     are stable for downward recursion and Neumann functions are stable
+     for upward recursion. It is
+     x2rev = 2.0/x,
+     J(N-1, x) = x2rev * N * J(N, x) - J(N+1, x),
+     Y(N+1, x) = x2rev * N * Y(N, x) - Y(N-1, x).  */
+
+  gfc_set_model_kind (x->ts.kind);
+
+  /* Get first recursion anchor.  */
+
+  mpfr_init (last1);
+  if (jn)
+    mpfr_jn (last1, n2, x->value.real, GFC_RND_MODE);
+  else
+    mpfr_yn (last1, n1, x->value.real, GFC_RND_MODE);
+
+  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_set (e->value.real, last1, GFC_RND_MODE);
+  if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
+    {
+      mpfr_clear (last1);
+      gfc_free_expr (e);
+      gfc_free_expr (result);
+      return &gfc_bad_expr;
+    }
+  gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
+
+  if (n1 == n2)
+    {
+      mpfr_clear (last1);
+      return result;
+    }
+ 
+  /* Get second recursion anchor.  */
+
+  mpfr_init (last2);
+  if (jn)
+    mpfr_jn (last2, n2-1, x->value.real, GFC_RND_MODE);
+  else
+    mpfr_yn (last2, n1+1, x->value.real, GFC_RND_MODE);
+
+  e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+  mpfr_set (e->value.real, last2, GFC_RND_MODE);
+  if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
+    {
+      mpfr_clear (last1);
+      mpfr_clear (last2);
+      gfc_free_expr (e);
+      gfc_free_expr (result);
+      return &gfc_bad_expr;
+    }
+  if (jn)
+     gfc_constructor_insert_expr (&result->value.constructor, e, &x->where, -2);
+  else 
+    gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
+
+  if (n1 + 1 == n2)
+    {
+      mpfr_clear (last1);
+      mpfr_clear (last2);
+      return result;
+    }
+
+  /* Start actual recursion.  */
+
+  mpfr_init (x2rev);
+  mpfr_ui_div (x2rev, 2.0, x->value.real, GFC_RND_MODE);
+ 
+  for (i = 2; i <= n2-n1; i++)
+    {
+      e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+      mpfr_mul_si (e->value.real, x2rev, jn ? (n2-i+1) : (n1+i-1),
+		   GFC_RND_MODE);
+      mpfr_mul (e->value.real, e->value.real, last2, GFC_RND_MODE);
+      mpfr_sub (e->value.real, e->value.real, last1, GFC_RND_MODE);
+
+      if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
+	goto error;
+
+      if (jn)
+	gfc_constructor_insert_expr (&result->value.constructor, e, &x->where,
+				     -i-1);
+      else
+	gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
+
+      mpfr_set (last1, last2, GFC_RND_MODE);
+      mpfr_set (last2, e->value.real, GFC_RND_MODE);
+    }
+
+  mpfr_clear (last1);
+  mpfr_clear (last2);
+  mpfr_clear (x2rev);
+  return result;
+
+error:
+  mpfr_clear (last1);
+  mpfr_clear (last2);
+  mpfr_clear (x2rev);
+  gfc_free_expr (e);
+  gfc_free_expr (result);
+  return &gfc_bad_expr;
+}
+
+
+gfc_expr *
+gfc_simplify_bessel_jn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
+{
+  return gfc_simplify_bessel_n2 (order1, order2, x, true);
+}
+
+
 gfc_expr *
 gfc_simplify_bessel_y0 (gfc_expr *x)
 {
@@ -1243,6 +1420,13 @@  gfc_simplify_bessel_yn (gfc_expr *order,
 }
 
 
+gfc_expr *
+gfc_simplify_bessel_yn2 (gfc_expr *order1, gfc_expr *order2, gfc_expr *x)
+{
+  return gfc_simplify_bessel_n2 (order1, order2, x, false);
+}
+
+
 gfc_expr *
 gfc_simplify_bit_size (gfc_expr *e)
 {
Index: gcc/testsuite/gfortran.dg/bessel_3.f90
===================================================================
--- gcc/testsuite/gfortran.dg/bessel_3.f90	(Revision 0)
+++ gcc/testsuite/gfortran.dg/bessel_3.f90	(Revision 0)
@@ -0,0 +1,18 @@ 
+! { dg-do compile }
+! { dg-options "-std=f2003 -Wimplicit-procedure" }
+!
+! PR fortran/36158 - Transformational BESSEL_JN/YN
+! PR fortran/33197 - F2008 math functions
+!
+IMPLICIT NONE
+print *, SIN (1.0)
+print *, BESSEL_J0(1.0) ! { dg-error "has no IMPLICIT type" })
+print *, BESSEL_J1(1.0) ! { dg-error "has no IMPLICIT type" }
+print *, BESSEL_JN(1,1.0) ! { dg-error "has no IMPLICIT type" }
+print *, BESSEL_JN(1,2,1.0) ! { dg-error "has no IMPLICIT type" }
+
+print *, BESSEL_Y0(1.0) ! { dg-error "has no IMPLICIT type" }
+print *, BESSEL_Y1(1.0) ! { dg-error "has no IMPLICIT type" }
+print *, BESSEL_YN(1,1.0) ! { dg-error "has no IMPLICIT type" }
+print *, BESSEL_YN(1,2,1.0) ! { dg-error "has no IMPLICIT type" }
+end
Index: gcc/testsuite/gfortran.dg/bessel_4.f90
===================================================================
--- gcc/testsuite/gfortran.dg/bessel_4.f90	(Revision 0)
+++ gcc/testsuite/gfortran.dg/bessel_4.f90	(Revision 0)
@@ -0,0 +1,24 @@ 
+! { dg-do compile }
+! { dg-options "-std=f2008" }
+!
+! PR fortran/36158 - Transformational BESSEL_JN/YN
+! PR fortran/33197 - F2008 math functions
+!
+implicit none
+! OK, elemental function:
+ print *, bessel_yn(1, [1.0, 2.0])
+ print *, bessel_yn([1, 2], 2.0)
+
+! Wrong, transformational function:
+! Does not pass check.c -- thus regarded as wrong generic function
+! and thus rejected with a slightly misleading error message
+ print *, bessel_yn(1, 2, [2.0, 3.0]) ! { dg-error "Too many arguments" }
+
+! Wrong in F2008: Negative argument, ok as GNU extension
+ print *, bessel_yn(-1, 3.0) ! { dg-error "Extension: Negative argument N " }
+
+! Wrong in F2008: Negative argument -- and no need for a GNU extension
+! Does not pass check.c -- thus regarded as wrong generic function
+! and thus rejected with a slightly misleading error message
+ print *, bessel_yn(-1, 2, 3.0) ! { dg-error "Too many arguments" }
+end
Index: gcc/testsuite/gfortran.dg/bessel_5.f90
===================================================================
--- gcc/testsuite/gfortran.dg/bessel_5.f90	(Revision 0)
+++ gcc/testsuite/gfortran.dg/bessel_5.f90	(Revision 0)
@@ -0,0 +1,86 @@ 
+! { dg-do run }
+! { dg-options "-Wall -fno-range-check" }
+!
+! PR fortran/36158 - Transformational BESSEL_JN/YN
+! PR fortran/33197 - F2008 math functions
+!
+! This is a dg-do run test as the middle end cannot simplify the
+! the scalarization of the elemental function (cf. PR 45305).
+!
+! -Wall has been specified to disabled -pedantic, which warns about the
+! negative order (GNU extension) to the order of the Bessel functions of
+! first and second kind.
+!
+
+implicit none
+integer :: i
+
+
+! Difference to mpfr_jn: 1 ULP
+
+if (any (abs (BESSEL_JN(2, 5, 2.457) - [(BESSEL_JN(i, 2.457), i = 2, 5)]) &
+          > epsilon(0.0))) then
+  print *, 'FAIL 1'
+  call abort()
+end if
+
+
+! Difference to mpfr_yn: 4 ULP
+
+if (any (abs (BESSEL_YN(2, 5, 2.457) - [(BESSEL_YN(i, 2.457), i = 2, 5)]) &
+         > epsilon(0.0)*4)) then
+  call abort()
+end if
+
+
+! Difference to mpfr_jn: 1 ULP
+
+if (any (abs (BESSEL_JN(0, 10, 4.457) &
+              - [ (BESSEL_JN(i, 4.457), i = 0, 10) ]) &
+         > epsilon(0.0))) then
+  call abort()
+end if
+
+
+! Difference to mpfr_yn: 192 ULP
+
+if (any (abs (BESSEL_YN(0, 10, 4.457) &
+              - [ (BESSEL_YN(i, 4.457), i = 0, 10) ]) &
+         > epsilon(0.0)*192)) then
+  call abort()
+end if
+
+
+! Difference to mpfr_jn: 0 ULP
+
+if (any (BESSEL_JN(0, 10, 0.0) /= [ (BESSEL_JN(i, 0.0), i = 0, 10) ])) &
+then
+  call abort()
+end if
+
+
+! Difference to mpfr_yn: 0 ULP
+
+if (any (BESSEL_YN(0, 10, 0.0) /= [ (BESSEL_YN(i, 0.0), i = 0, 10) ])) &
+then
+  call abort()
+end if
+
+
+! Difference to mpfr_jn: 1 ULP
+
+if (any (abs (BESSEL_JN(0, 10, 1.0) &
+              - [ (BESSEL_JN(i, 1.0), i = 0, 10) ]) &
+         > epsilon(0.0)*1)) then
+ call abort()
+end if
+
+! Difference to mpfr_yn: 32 ULP
+
+if (any (abs (BESSEL_YN(0, 10, 1.0) &
+              - [ (BESSEL_YN(i, 1.0), i = 0, 10) ]) &
+         > epsilon(0.0)*32)) then
+  call abort()
+end if
+
+end