From patchwork Tue Aug 17 12:37:54 2010 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Tobias Burnus X-Patchwork-Id: 61890 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 A319FB70AE for ; Tue, 17 Aug 2010 22:38:17 +1000 (EST) Received: (qmail 30363 invoked by alias); 17 Aug 2010 12:38:11 -0000 Received: (qmail 30101 invoked by uid 22791); 17 Aug 2010 12:38:06 -0000 X-SWARE-Spam-Status: No, hits=-0.9 required=5.0 tests=AWL, BAYES_20, RCVD_IN_DNSWL_NONE 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 12:37:59 +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 DAF803D813; Tue, 17 Aug 2010 14:37:54 +0200 (CEST) Message-ID: <4C6A82A2.6020909@net-b.de> Date: Tue, 17 Aug 2010 14:37:54 +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: gcc patches , gfortran Subject: [Patch, Fortran] PR36158 - Add transformational version of BESSEL_JN intrinsic 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 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. 2010-08-17 Tobias Burnus 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 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