diff mbox

[Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic

Message ID 4C6A82A2.6020909@net-b.de
State New
Headers show

Commit Message

Tobias Burnus Aug. 17, 2010, 12:37 p.m. UTC 
Fortran 2008 has
   BESSEL_JN(N, X)
as elemental function (which is implemented) and
   BESSEL_JN(N1, N2, X)
as transformational function - which this patch adds.

Note: This patch only does the compile-time implementation - a run-time 
implementation is still pending.

Build and regtested on x86-64-linux.
OK for the trunk? (Pending the issue below.)

Note: I get the following failures, which seem to be unrelated:
a) gfortran.dg/widechar_intrinsics_5.f90: Fails (coredump) at -O1 but 
not at -O0
b) gfortran.dg/trim_optimize_1.f90: Fails in the dump for having more 
than 2 memmove
On gcc-testresults, I do not see those failures, but the builds there 
are also a tad older. I am currently bootstrapping to make sure. (b) 
seems to be a FE issue (dependency.c, frontend-passes.c) - but I do not 
have any modifications there; (a) looks like a FE->ME issue, but I do 
not see how this patch could affect those.

Tobias

PS: The array_memcpy_3.f90 issue (PR45266) has been solved by Richard by 
adjusting the pattern.

Comments

Daniel Kraft Aug. 17, 2010, 1:23 p.m. UTC | #1 
Tobias Burnus wrote:
>  Fortran 2008 has
>   BESSEL_JN(N, X)
> as elemental function (which is implemented) and
>   BESSEL_JN(N1, N2, X)
> as transformational function - which this patch adds.
> 
> Note: This patch only does the compile-time implementation - a run-time 
> implementation is still pending.
> 
> Build and regtested on x86-64-linux.
> OK for the trunk? (Pending the issue below.)

Ok.  But:

+      if (n < 0 && gfc_notify_std (GFC_STD_GNU, "Extension: Negativ 
argument "
+				    "N at %L", &order->where) == FAILURE)

It's "negative" in English, I think.  The same is later, also.

> Note: I get the following failures, which seem to be unrelated:
> a) gfortran.dg/widechar_intrinsics_5.f90: Fails (coredump) at -O1 but 
> not at -O0
> b) gfortran.dg/trim_optimize_1.f90: Fails in the dump for having more 
> than 2 memmove
> On gcc-testresults, I do not see those failures, but the builds there 
> are also a tad older. I am currently bootstrapping to make sure. (b) 
> seems to be a FE issue (dependency.c, frontend-passes.c) - but I do not 
> have any modifications there; (a) looks like a FE->ME issue, but I do 
> not see how this patch could affect those.

I see (b) but not (a) on my system with the build of my last ASSOCIATE 
patch, but there has been some partial svn update so it may be wrong.

Yours,
Daniel
Steve Kargl Aug. 17, 2010, 2:14 p.m. UTC | #2 
On Tue, Aug 17, 2010 at 03:23:14PM +0200, Daniel Kraft wrote:
> Tobias Burnus wrote:
> > Fortran 2008 has
> >  BESSEL_JN(N, X)
> >as elemental function (which is implemented) and
> >  BESSEL_JN(N1, N2, X)
> >as transformational function - which this patch adds.
> >
> >Note: This patch only does the compile-time implementation - a run-time 
> >implementation is still pending.
> >
> >Build and regtested on x86-64-linux.
> >OK for the trunk? (Pending the issue below.)
> 
> Ok.  But:
> 

Not OK. :-)

In check.c, 

+  if (scalar_check (n2, 0) == FAILURE)
+    return FAILURE;

The 0 should 1.

+
+  if (type_check (x, 2, BT_REAL) == FAILURE)
+    return FAILURE;
+  if (scalar_check (x, 0) == FAILURE)
+    return FAILURE;

The 0 should be 2.

Also note that if n1 and n2 are to be nonnegative.
So, one should add,

   if (nonnegative_check ('n1', n1) == FAILURE)
     return FAILURE;

   if (nonnegative_check ('n2', n2) == FAILURE)
     return FAILURE;


In the simplification route

+  if (order->expr_type == EXPR_CONSTANT)
+    {
+      n = mpz_get_si (order->value.integer);
+      if (n < 0 && gfc_notify_std (GFC_STD_GNU, "Extension: Negativ argument "
+                                   "N at %L", &order->where) == FAILURE)
+       return &gfc_bad_expr;
+    }

   if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
     return NULL;

Please move the 2nd if-statement to the top of the function.  There's
no reason to possibly issue an error for n < 0, if once that is fixed
simplification terminates due to a non-constant x.

Finally, note that the mpfr_{jn,yn} routines are called n2 - n1 + 1
times, which may be very slow for large n2 - n1 + 1.  Bessel functions
are stable for downward recursion, and Neumann functions are stable
for upward recursion.  It may be better to use recursion, which 
involves two evaluations of mfpr_{jn,yn} and then some simple
multiplications and additions.
Tobias Burnus Aug. 17, 2010, 2:50 p.m. UTC | #3 
On 08/17/2010 04:14 PM, Steve Kargl wrote:
> Not OK. :-)
> In check.c,
> +  if (scalar_check (n2, 0) == FAILURE)
> +    return FAILURE;
> The 0 should 1.
Thanks for spotting this.

> Also note that if n1 and n2 are to be nonnegative.
> So, one should add,
>
>     if (nonnegative_check ('n1', n1) == FAILURE)
>       return FAILURE;
>     if (nonnegative_check ('n2', n2) == FAILURE)
>       return FAILURE;

I disagree. Negative values are perfectly valid thus I would prefer 
adding such a restriction only for -std=f2003.

(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.)

> In the simplification route
>
> +  if (order->expr_type == EXPR_CONSTANT)
> +    {
> +      n = mpz_get_si (order->value.integer);
> +      if (n<  0&&  gfc_notify_std (GFC_STD_GNU, "Extension: Negativ argument "
> +                                   "N at %L",&order->where) == FAILURE)
> +       return&gfc_bad_expr;
> +    }
>
>     if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
>       return NULL;
>
> Please move the 2nd if-statement to the top of the function.  There's
> no reason to possibly issue an error for n<  0, if once that is fixed
> simplification terminates due to a non-constant x.

Well, the error is useful, unless one already puts the error into 
check.c. But I agree that one can move the check into check.c.

> Finally, note that the mpfr_{jn,yn} routines are called n2 - n1 + 1
> times, which may be very slow for large n2 - n1 + 1.  Bessel functions
> are stable for downward recursion, and Neumann functions are stable
> for upward recursion.  It may be better to use recursion, which
> involves two evaluations of mfpr_{jn,yn} and then some simple
> multiplications and additions

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)

Tobias
Steve Kargl Aug. 17, 2010, 3:36 p.m. UTC | #4 
On Tue, Aug 17, 2010 at 04:50:09PM +0200, Tobias Burnus wrote:
>  On 08/17/2010 04:14 PM, Steve Kargl wrote:
> >Not OK. :-)
> >In check.c,
> >+  if (scalar_check (n2, 0) == FAILURE)
> >+    return FAILURE;
> >The 0 should 1.
> Thanks for spotting this.
> 
> >Also note that if n1 and n2 are to be nonnegative.
> >So, one should add,
> >
> >    if (nonnegative_check ('n1', n1) == FAILURE)
> >      return FAILURE;
> >    if (nonnegative_check ('n2', n2) == FAILURE)
> >      return FAILURE;
> 
> 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.

+  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. 

> (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.

> >In the simplification route
> >
> >+  if (order->expr_type == EXPR_CONSTANT)
> >+    {
> >+      n = mpz_get_si (order->value.integer);
> >+      if (n<  0&&  gfc_notify_std (GFC_STD_GNU, "Extension: Negativ 
> >argument "
> >+                                   "N at %L",&order->where) == FAILURE)
> >+       return&gfc_bad_expr;
> >+    }
> >
> >    if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
> >      return NULL;
> >
> >Please move the 2nd if-statement to the top of the function.  There's
> >no reason to possibly issue an error for n<  0, if once that is fixed
> >simplification terminates due to a non-constant x.
> 
> Well, the error is useful, unless one already puts the error into 
> check.c. But I agree that one can move the check into check.c.

Yes, the check belong in check.c.

> >Finally, note that the mpfr_{jn,yn} routines are called n2 - n1 + 1
> >times, which may be very slow for large n2 - n1 + 1.  Bessel functions
> >are stable for downward recursion, and Neumann functions are stable
> >for upward recursion.  It may be better to use recursion, which
> >involves two evaluations of mfpr_{jn,yn} and then some simple
> >multiplications and additions
> 
> 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 you have Abromowitz and Stegun, look up Miller's
algorithm for computing Bessel functions.  The new NIST
rewrite has

http://dlmf.nist.gov/10.74#iv

The guts of my bessel function routine has

   if (x == 0.e0_knd) then  ! Avoid division by zero
      j = 0
      j(0) = 1
   else
      !  
      ! Do the downward recursion.
      !
      jmx1 = besjn(n + 2, x)
      jmx  = besjn(n + 1, x)
      ix = 2 / x
      do i = n + 1, 1, -1
         j(i - 1) = i * ix * jmx - jmx1
         jmx1 = jmx
         jmx = j(i - 1)
      end do

If the initial values of J(N,x) and J(N+1, x) are accurate
approximations then the above generates the required sequence.
If the values are inaccurate, the generated sequence is 
porportional to the desired sequence of Bessel functions.
One then to normalize the sequence, for my routine I do

      jmx = cj0(x)   ! J(0)
      jmx1 = cj1(x)  ! J(1)
      if (abs(jmx) >= abs(jmx1)) then
         jmx = jmx  / j(0)
      else
         jmx = jmx1 / j(1)
      end if
      j = jmx * j
   end if
diff mbox

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/check.c
===================================================================
--- gcc/fortran/check.c	(revision 163302)
+++ gcc/fortran/check.c	(working copy)
@@ -892,6 +892,31 @@  gfc_check_besn (gfc_expr *n, gfc_expr *x
 }
 
 
+/* 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 (type_check (n2, 1, BT_INTEGER) == FAILURE)
+    return FAILURE;
+  if (scalar_check (n2, 0) == FAILURE)
+    return FAILURE;
+
+  if (type_check (x, 2, BT_REAL) == FAILURE)
+    return FAILURE;
+  if (scalar_check (x, 0) == FAILURE)
+    return FAILURE;
+
+  return SUCCESS;
+}
+
+
 gfc_try
 gfc_check_bitfcn (gfc_expr *i, gfc_expr *pos)
 {
Index: gcc/fortran/gfortran.h
===================================================================
--- gcc/fortran/gfortran.h	(revision 163302)
+++ gcc/fortran/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: gcc/fortran/intrinsic.c
===================================================================
--- gcc/fortran/intrinsic.c	(revision 163302)
+++ gcc/fortran/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: gcc/fortran/intrinsic.h
===================================================================
--- gcc/fortran/intrinsic.h	(revision 163302)
+++ gcc/fortran/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: gcc/fortran/intrinsic.texi
===================================================================
--- gcc/fortran/intrinsic.texi	(revision 163302)
+++ gcc/fortran/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: gcc/fortran/simplify.c
===================================================================
--- gcc/fortran/simplify.c	(revision 163302)
+++ gcc/fortran/simplify.c	(working copy)
@@ -1183,12 +1183,19 @@  gfc_expr *
 gfc_simplify_bessel_jn (gfc_expr *order, gfc_expr *x)
 {
   gfc_expr *result;
-  long n;
+  long n = 0;
+
+  if (order->expr_type == EXPR_CONSTANT)
+    {
+      n = mpz_get_si (order->value.integer);
+      if (n < 0 && gfc_notify_std (GFC_STD_GNU, "Extension: Negativ argument "
+				    "N at %L", &order->where) == FAILURE)
+	return &gfc_bad_expr;
+    }
 
   if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  n = mpz_get_si (order->value.integer);
   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
   mpfr_jn (result->value.real, n, x->value.real, GFC_RND_MODE);
 
@@ -1196,6 +1203,72 @@  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;
+  long n1 = 0, n2 = 0;
+  int i;
+
+  if (order1->expr_type == EXPR_CONSTANT)
+    {
+      n1 = mpz_get_si (order1->value.integer);
+      if (n1 < 0 && gfc_notify_std (GFC_STD_GNU, "Extension: Negativ argument "
+				    "N1 at %L", &order1->where) == FAILURE)
+	return &gfc_bad_expr;
+    }
+
+  if (order2->expr_type == EXPR_CONSTANT)
+    {
+      n2 = mpz_get_si (order2->value.integer);
+      if (n2 < 0 && gfc_notify_std (GFC_STD_GNU, "Extension: Negativ argument "
+				    "N2 at %L", &order2->where) == FAILURE)
+	return &gfc_bad_expr;
+    }
+
+  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;
+    }
+
+  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));
+
+  for (i = n1; i <= n2; ++i)
+    {
+      gfc_expr *e;
+     
+      e = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
+
+      if (jn)
+	mpfr_jn (e->value.real, i, x->value.real, GFC_RND_MODE);
+      else
+	mpfr_yn (e->value.real, i, x->value.real, GFC_RND_MODE);
+
+      if (range_check (e, jn ? "BESSEL_JN" : "BESSEL_YN") == &gfc_bad_expr)
+	return &gfc_bad_expr;
+      gfc_constructor_append_expr (&result->value.constructor, e, &x->where);
+    }
+
+  return 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)
 {
@@ -1230,12 +1303,19 @@  gfc_expr *
 gfc_simplify_bessel_yn (gfc_expr *order, gfc_expr *x)
 {
   gfc_expr *result;
-  long n;
+  long n = 0;
+
+  if (order->expr_type == EXPR_CONSTANT)
+    {
+      n = mpz_get_si (order->value.integer);
+      if (n < 0 && gfc_notify_std (GFC_STD_GNU, "Extension: Negativ argument "
+				    "N at %L", &order->where) == FAILURE)
+	return &gfc_bad_expr;
+    }
 
   if (x->expr_type != EXPR_CONSTANT || order->expr_type != EXPR_CONSTANT)
     return NULL;
 
-  n = mpz_get_si (order->value.integer);
   result = gfc_get_constant_expr (x->ts.type, x->ts.kind, &x->where);
   mpfr_yn (result->value.real, n, x->value.real, GFC_RND_MODE);
 
@@ -1243,6 +1323,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,16 @@ 
+! { 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:
+ 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: Negativ argument N " }
+ print *, bessel_yn(-1, 2, 3.0) ! { dg-error "Extension: Negativ argument N1" }
+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,35 @@ 
+! { dg-do run }
+! { dg-options "-Wall" }
+!
+! 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
+if (any (abs (BESSEL_JN(2, 5, 2.457) - BESSEL_JN([2,3,4,5], 2.457)) &
+          > epsilon(0.0))) &
+  call abort()
+if (any (abs (BESSEL_YN(2, 5, 2.457) - BESSEL_YN([2,3,4,5], 2.457)) &
+         > epsilon(0.0))) &
+  call abort()
+if (any (abs (BESSEL_JN(-1, 1, 2.457) &
+              - [ BESSEL_JN(-1, 2.457), BESSEL_JN(0, 2.457), BESSEL_JN(1, 2.457) ]) &
+         > epsilon(0.0))) &
+  call abort()
+if (any (abs (BESSEL_YN(-1, 1, 2.457) &
+              - [ BESSEL_YN(-1, 2.457), BESSEL_YN(0, 2.457), BESSEL_YN(1, 2.457) ]) &
+         > epsilon(0.0))) &
+  call abort()
+if (any (abs (BESSEL_JN(0, 1, 2.457) - [BESSEL_J0(2.457), BESSEL_J1(2.457)]) &
+    > epsilon(0.0))) &
+  call abort()
+if (any (abs (BESSEL_YN(0, 1, 2.457) - [BESSEL_Y0(2.457), BESSEL_Y1(2.457)]) &
+         > epsilon(0.0))) &
+  call abort()
+end