OpenE2K
/
gcc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

re PR fortran/36158 (Transformational function BESSEL_YN(n1,n2,x) and BESSEL_JN missing) 

2010-08-19  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-19  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.

From-SVN: r163364
This commit is contained in:
Tobias Burnus 2010-08-19 09:28:17 +02:00 committed by Tobias Burnus
parent 771c5727a0
commit 29698e0f2f
11 changed files with 425 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
gcc/fortran/ChangeLog
View File

@ -1,3 +1,20 @@
2010-08-19 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-19 Jerry DeLisle <jvdelisle@gcc.gnu.org>
PR fortran/41859

36
gcc/fortran/check.c
View File

@ -884,6 +884,14 @@ 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;
@ -892,6 +900,34 @@ 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 (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;
}
gfc_try
gfc_check_bitfcn (gfc_expr *i, gfc_expr *pos)
{

4
gcc/fortran/gfortran.h
View File

@ -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;

10
gcc/fortran/intrinsic.c
View File

@ -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,

3
gcc/fortran/intrinsic.h
View File

@ -40,6 +40,7 @@ gfc_try gfc_check_associated (gfc_expr *, 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_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 *);

46
gcc/fortran/intrinsic.texi
View File

@ -1630,29 +1630,41 @@ 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} is 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{X} @tab Shall be a scalar or an array of type @code{REAL}.
@item @var{N1} @tab Shall be a non-negative scalar of type @code{INTEGER}.
@item @var{N2} @tab Shall be a non-negative 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}:
The return value is a scalar of type @code{REAL}. It has the same
kind as @var{X}.
@item @emph{Note}:
The transformational function uses a recurrance algorithm which might,
for some values of @var{X}, lead to different results than calls to
the elemental function.
@item @emph{Example}:
@smallexample
program test_besjn
@ -1778,29 +1790,41 @@ 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} is 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 non-negative scalar of type @code{INTEGER}.
@item @var{N2} @tab Shall be a non-negative scalar of type @code{INTEGER}.
@item @var{X} @tab Shall be a scalar or an array of type @code{REAL};
for @code{BESSEL_YN(N1, N2, X)} it shall be scalar.
@end multitable
@item @emph{Return value}:
The return value is a scalar of type @code{REAL}. It has the same
kind as @var{X}.
@item @emph{Note}:
The transformational function uses a recurrance algorithm which might,
for some values of @var{X}, lead to different results than calls to
the elemental function.
@item @emph{Example}:
@smallexample
program test_besyn

185
gcc/fortran/simplify.c
View File

@ -1196,6 +1196,184 @@ gfc_simplify_bessel_jn (gfc_expr *order, gfc_expr *x)
}
/* 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).
Cf. http://dlmf.nist.gov/10.74#iv and http://dlmf.nist.gov/10.6#E1 */
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, 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 +1421,13 @@ gfc_simplify_bessel_yn (gfc_expr *order, gfc_expr *x)
}
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)
{

8
gcc/testsuite/ChangeLog
View File

@ -1,3 +1,11 @@
2010-08-19 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.
2010-08-19 Janus Weil <janus@gcc.gnu.org>
PR fortran/45290

18
gcc/testsuite/gfortran.dg/bessel_3.f90 Normal file
View File

@ -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

24
gcc/testsuite/gfortran.dg/bessel_4.f90 Normal file
View File

@ -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

86
gcc/testsuite/gfortran.dg/bessel_5.f90 Normal file
View File

@ -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 epsilon
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 epsilon
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 epsilon
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 epsilon
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: None. (Special case: X = 0.0)
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: None. (Special case: X = 0.0)
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 epsilon
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 epsilon
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