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

Message ID 4C6ABFE0.7030908@net-b.de New show

## Commit Message

Tobias Burnus Aug. 17, 2010, 4:59 p.m.
```On 08/17/2010 05:36 PM, Steve Kargl wrote:
>> I disagree. Negative values are perfectly valid thus I would prefer
>> adding such a restriction only for -std=f2003.
> Then you may need to fix your simplicification function.

Let's do it differently: For the existing elemental function, we allow
non-negative values but for the transformational one, we don't.

> +  mpz_init_set_ui (result->shape[0], MAX (n2-n1+1, 0));
> +
> +  for (i = n1; i<= n2; ++i)
> +    {
>
> What happens if n2 = -4 and n1 = 1?  Note, the standard
> does not specify that n2>  n1.  It is implied by
>
>      Case (ii): The result of BESSEL JN (N1, N2, X) is a rank-one
>        array with extent MAX (N2-N1+1, 0).
>
> but there is no requirement for this ordering.  It may be prudent to
> add a check that n2>  n1.

I think that's wrong. First, n2 == n1 is perfectly valid; secondly, if
n2 < n1, one gets a zero-sized array, which might be not terribly
useful, but is perfectly valid.

>> (J(-m,x) = (-1)**m * J(m,x) -- and analogously for YN - at least for
>> integral m; using the Gamma function for negative m won't work, as it
>> becomes +/-infinite for negative integer values). Thus, there is no need
>> for negative "m", but also no need to reject them. As gfortran allowed
>> negative "m" so far, I think it should continue to do so for -std=gnu.)
> I would prefer to remove this extension.
Why?

>> Do you mean the following algorithm?
>> 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
> Yes.
>     if (x == 0.e0_knd) then  ! Avoid division by zero
>        j = 0
The x == 0 is a good point!

How about the attached patch? The patch seems to work correctly for JN
and x == 0  and for YN (all x); the result agrees with the GLIBC result
up to ULP.

The problem I still have is with JN. My feeling is that
gfc_constructor_insert_expr does not correctly work. For
print *, BESSEL_JN(0, 1, 1.0)
print *, (BESSEL_JN(i, 1.0), i =0, 1)
I get:
0.76519769       0.0000000
0.76519769      0.44005057
using
gfc_constructor_insert_expr (&result->value.constructor, e,
&x->where, 0);
if I change the "0" to "1", I get:
0.44005057      0.76519769
thus the value is correctly set - but somehow inserting the constructor
as first element does not work.

Any idea?

Tobias

PS: Steve, are you interested in implementing a TREE or libgfortran
version of this algorithm? I think using a libgfortran version makes
more sense than implementing it as TREE; what do you think?

PPS: See first email for some test cases or see below.

print *, BESSEL_JN(0, 5, 0.0)
print *, (BESSEL_JN(i, 0.0), i =0, 5)
print *,
print *, BESSEL_YN(0, 5, 0.0)
print *, (BESSEL_YN(i, 0.0), i =0, 5)
print *,
print *, BESSEL_JN(0, 1, 1.0)
print *, (BESSEL_JN(i, 1.0), i =0, 1)
print *, ''
print *, BESSEL_YN(0, 4 2.5)
print *, (BESSEL_YN(i, 2.5), i =0, 4)
end
```

Steve Kargl Aug. 17, 2010, 7:01 p.m. | #1
```On Tue, Aug 17, 2010 at 06:59:12PM +0200, Tobias Burnus wrote:
>  On 08/17/2010 05:36 PM, Steve Kargl wrote:
> >>I disagree. Negative values are perfectly valid thus I would prefer
> >>adding such a restriction only for -std=f2003.
> >Then you may need to fix your simplicification function.
>
> Let's do it differently: For the existing elemental function, we allow
> non-negative values but for the transformational one, we don't.
>
> >+  mpz_init_set_ui (result->shape[0], MAX (n2-n1+1, 0));
> >+
> >+  for (i = n1; i<= n2; ++i)
> >+    {
> >
> >What happens if n2 = -4 and n1 = 1?  Note, the standard
> >does not specify that n2>  n1.  It is implied by
> >
> >     Case (ii): The result of BESSEL JN (N1, N2, X) is a rank-one
> >       array with extent MAX (N2-N1+1, 0).
> >
> >but there is no requirement for this ordering.  It may be prudent to
> >add a check that n2>  n1.
>
> I think that's wrong. First, n2 == n1 is perfectly valid; secondly, if
> n2 < n1, one gets a zero-sized array, which might be not terribly
> useful, but is perfectly valid.

Whoops.  I misread the condition in the for-loop,
and thought that you would enter an infinite loop
stomping on memory.

> >>(J(-m,x) = (-1)**m * J(m,x) -- and analogously for YN - at least for
> >>integral m; using the Gamma function for negative m won't work, as it
> >>becomes +/-infinite for negative integer values). Thus, there is no need
> >>for negative "m", but also no need to reject them. As gfortran allowed
> >>negative "m" so far, I think it should continue to do so for -std=gnu.)
> >I would prefer to remove this extension.
> Why?

Experience.  I've been computing/using J() and Y() for more
than 23 years.  Getting a correct and accurate value is
is difficult.  The functions in libm have horrible accuracy
near their zeros.

> >>Do you mean the following algorithm?
> >>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
> >Yes.
> >    if (x == 0.e0_knd) then  ! Avoid division by zero
> >       j = 0
> The x == 0 is a good point!
>
> How about the attached patch?

I'll look over the patch in more detail in a few hours.

> PS: Steve, are you interested in implementing a TREE or libgfortran
> version of this algorithm? I think using a libgfortran version makes
> more sense than implementing it as TREE; what do you think?

I think we want to use TREE only because mpfr_{jn,yn} will
give accurate values for the two initial values in the
recursion scheme.  As noted above, the libm routines give
horrible accuracy near their zeros.  For j0(), I found
over a million ULP at its 1st zero.  If you're interested
in the details

http://www.freebsd.org/cgi/query-pr.cgi?pr=standards/142803
```

## Patch

```Index: intrinsic.c
===================================================================
--- intrinsic.c	(revision 163305)
+++ intrinsic.c	(working copy)
@@ -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: intrinsic.h
===================================================================
--- intrinsic.h	(revision 163305)
+++ intrinsic.h	(working copy)
@@ -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: decl.c
===================================================================
--- decl.c	(revision 163305)
+++ decl.c	(working copy)
@@ -1312,9 +1312,10 @@  add_init_expr_to_sym (const char *name,
}

/* Check if the assignment can happen. This has to be put off
-	 until later for a derived type variable.  */
+	 until later for derived type variables and procedure pointers.  */
if (sym->ts.type != BT_DERIVED && init->ts.type != BT_DERIVED
&& sym->ts.type != BT_CLASS && init->ts.type != BT_CLASS
+	  && !sym->attr.proc_pointer
&& gfc_check_assign_symbol (sym, init) == FAILURE)
return FAILURE;

@@ -1652,6 +1653,48 @@  gfc_match_null (gfc_expr **result)
}

+/* Match the initialization expr for a data pointer or procedure pointer.  */
+
+static match
+match_pointer_init (gfc_expr **init, int procptr)
+{
+  match m;
+
+  if (gfc_pure (NULL) && gfc_state_stack->state != COMP_DERIVED)
+    {
+      gfc_error ("Initialization of pointer at %C is not allowed in "
+		 "a PURE procedure");
+      return MATCH_ERROR;
+    }
+
+  /* Match NULL() initilization.  */
+  m = gfc_match_null (init);
+  if (m != MATCH_NO)
+    return m;
+
+  /* Match non-NULL initialization.  */
+  gfc_matching_procptr_assignment = procptr;
+  m = gfc_match_rvalue (init);
+  gfc_matching_procptr_assignment = 0;
+  if (m == MATCH_ERROR)
+    return MATCH_ERROR;
+  else if (m == MATCH_NO)
+    {
+      gfc_error ("Error in pointer initialization at %C");
+      return MATCH_ERROR;
+    }
+
+  if (!procptr)
+    gfc_resolve_expr (*init);
+
+  if (gfc_notify_std (GFC_STD_F2008, "Fortran 2008: non-NULL pointer "
+		      "initialization at %C") == FAILURE)
+    return MATCH_ERROR;
+
+  return MATCH_YES;
+}
+
+
/* Match a variable name with an optional initializer.  When this
subroutine is called, a variable is expected to be parsed next.
Depending on what is happening at the moment, updates either the
@@ -1899,23 +1942,9 @@  variable_decl (int elem)
goto cleanup;
}

-	  m = gfc_match_null (&initializer);
-	  if (m == MATCH_NO)
-	    {
-	      gfc_error ("Pointer initialization requires a NULL() at %C");
-	      m = MATCH_ERROR;
-	    }
-
-	  if (gfc_pure (NULL) && gfc_state_stack->state != COMP_DERIVED)
-	    {
-	      gfc_error ("Initialization of pointer at %C is not allowed in "
-			 "a PURE procedure");
-	      m = MATCH_ERROR;
-	    }
-
+	  m = match_pointer_init (&initializer, 0);
if (m != MATCH_YES)
goto cleanup;
-
}
else if (gfc_match_char ('=') == MATCH_YES)
{
@@ -4675,20 +4704,7 @@  match_procedure_decl (void)
goto cleanup;
}

-	  m = gfc_match_null (&initializer);
-	  if (m == MATCH_NO)
-	    {
-	      gfc_error ("Pointer initialization requires a NULL() at %C");
-	      m = MATCH_ERROR;
-	    }
-
-	  if (gfc_pure (NULL))
-	    {
-	      gfc_error ("Initialization of pointer at %C is not allowed in "
-			 "a PURE procedure");
-	      m = MATCH_ERROR;
-	    }
-
+	  m = match_pointer_init (&initializer, 1);
if (m != MATCH_YES)
goto cleanup;

@@ -4815,18 +4831,7 @@  match_ppc_decl (void)

if (gfc_match (" =>") == MATCH_YES)
{
-	  m = gfc_match_null (&initializer);
-	  if (m == MATCH_NO)
-	    {
-	      gfc_error ("Pointer initialization requires a NULL() at %C");
-	      m = MATCH_ERROR;
-	    }
-	  if (gfc_pure (NULL))
-	    {
-	      gfc_error ("Initialization of pointer at %C is not allowed in "
-			 "a PURE procedure");
-	      m = MATCH_ERROR;
-	    }
+	  m = match_pointer_init (&initializer, 1);
if (m != MATCH_YES)
{
gfc_free_expr (initializer);
Index: gfortran.h
===================================================================
--- gfortran.h	(revision 163305)
+++ gfortran.h	(working copy)
@@ -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: check.c
===================================================================
--- check.c	(revision 163305)
+++ check.c	(working copy)
@@ -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: intrinsic.texi
===================================================================
--- intrinsic.texi	(revision 163305)
+++ intrinsic.texi	(working copy)
@@ -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: simplify.c
===================================================================
--- simplify.c	(revision 163305)
+++ simplify.c	(working copy)
@@ -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);
+    }
+
+  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) == 1, JN(N > 1, 0.0) == 0; and
+     YN(N, 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);
+	}
+
+      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);
+    }
+  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);
+    }
+  if (jn)
+    gfc_constructor_insert_expr (&result->value.constructor, e, &x->where, 0);
+  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,
+				     0);
+      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);
+}
+
+
+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)
{

```