From patchwork Tue Aug 17 16:59:12 2010 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Tobias Burnus X-Patchwork-Id: 61937 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) by ozlabs.org (Postfix) with SMTP id DEA0DB6F14 for ; Wed, 18 Aug 2010 02:59:34 +1000 (EST) Received: (qmail 23379 invoked by alias); 17 Aug 2010 16:59:29 -0000 Received: (qmail 23353 invoked by uid 22791); 17 Aug 2010 16:59:26 -0000 X-SWARE-Spam-Status: No, hits=-1.1 required=5.0 tests=AWL, BAYES_00, KAM_STOCKGEN, RCVD_IN_DNSWL_NONE, TW_BG X-Spam-Check-By: sourceware.org Received: from mx01.qsc.de (HELO mx01.qsc.de) (213.148.129.14) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 17 Aug 2010 16:59:19 +0000 Received: from [192.168.178.22] (port-92-204-42-5.dynamic.qsc.de [92.204.42.5]) by mx01.qsc.de (Postfix) with ESMTP id 912033D9B6; Tue, 17 Aug 2010 18:59:12 +0200 (CEST) Message-ID: <4C6ABFE0.7030908@net-b.de> Date: Tue, 17 Aug 2010 18:59:12 +0200 From: Tobias Burnus User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv:1.9.2.7) Gecko/20100714 SUSE/3.1.1 Thunderbird/3.1.1 MIME-Version: 1.0 To: Steve Kargl CC: Daniel Kraft , gcc patches , gfortran Subject: Re: [Patch, Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic References: <4C6A82A2.6020909@net-b.de> <4C6A8D42.2000609@domob.eu> <20100817141438.GA35172@troutmask.apl.washington.edu> <4C6AA1A1.2040303@net-b.de> <20100817153611.GA35578@troutmask.apl.washington.edu> In-Reply-To: <20100817153611.GA35578@troutmask.apl.washington.edu> Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Archive: List-Post: List-Help: Sender: gcc-patches-owner@gcc.gnu.org Delivered-To: mailing list gcc-patches@gcc.gnu.org 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 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); + 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) == 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); + 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, 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); + 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) {