From: Tobias Burnus Date: Thu, 19 Aug 2010 07:28:17 +0000 (+0200) Subject: re PR fortran/36158 (Transformational function BESSEL_YN(n1,n2,x) and BESSEL_JN missing) X-Git-Tag: releases/gcc-4.6.0~4941 X-Git-Url: http://git.ipfire.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=29698e0f2fa1c77242782bc5c8c6a9327b8d32cf;p=thirdparty%2Fgcc.git re PR fortran/36158 (Transformational function BESSEL_YN(n1,n2,x) and BESSEL_JN missing) 2010-08-19 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-19 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. From-SVN: r163364 --- diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog index 24e20843acbe..2c6e6f65a775 100644 --- a/gcc/fortran/ChangeLog +++ b/gcc/fortran/ChangeLog @@ -1,3 +1,20 @@ +2010-08-19 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-19 Jerry DeLisle PR fortran/41859 diff --git a/gcc/fortran/check.c b/gcc/fortran/check.c index ad040f1d2cc1..36efffa6dfa3 100644 --- a/gcc/fortran/check.c +++ b/gcc/fortran/check.c @@ -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) { diff --git a/gcc/fortran/gfortran.h b/gcc/fortran/gfortran.h index 89a8e504711c..5ca248801797 100644 --- a/gcc/fortran/gfortran.h +++ b/gcc/fortran/gfortran.h @@ -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; diff --git a/gcc/fortran/intrinsic.c b/gcc/fortran/intrinsic.c index cdfcca7bcf89..3751167e75b7 100644 --- a/gcc/fortran/intrinsic.c +++ b/gcc/fortran/intrinsic.c @@ -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, diff --git a/gcc/fortran/intrinsic.h b/gcc/fortran/intrinsic.h index 23272a8cd0dc..7780ebca232f 100644 --- a/gcc/fortran/intrinsic.h +++ b/gcc/fortran/intrinsic.h @@ -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 *); diff --git a/gcc/fortran/intrinsic.texi b/gcc/fortran/intrinsic.texi index 566050062a1c..3c82ffc2c142 100644 --- a/gcc/fortran/intrinsic.texi +++ b/gcc/fortran/intrinsic.texi @@ -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 diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c index b47f8ccc393d..d1e94af4db2b 100644 --- a/gcc/fortran/simplify.c +++ b/gcc/fortran/simplify.c @@ -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) { diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog index d033f9a9d57a..db26a97676ec 100644 --- a/gcc/testsuite/ChangeLog +++ b/gcc/testsuite/ChangeLog @@ -1,3 +1,11 @@ +2010-08-19 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. + 2010-08-19 Janus Weil PR fortran/45290 diff --git a/gcc/testsuite/gfortran.dg/bessel_3.f90 b/gcc/testsuite/gfortran.dg/bessel_3.f90 new file mode 100644 index 000000000000..271768dd890f --- /dev/null +++ b/gcc/testsuite/gfortran.dg/bessel_3.f90 @@ -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 diff --git a/gcc/testsuite/gfortran.dg/bessel_4.f90 b/gcc/testsuite/gfortran.dg/bessel_4.f90 new file mode 100644 index 000000000000..7da1bf9aa17d --- /dev/null +++ b/gcc/testsuite/gfortran.dg/bessel_4.f90 @@ -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 diff --git a/gcc/testsuite/gfortran.dg/bessel_5.f90 b/gcc/testsuite/gfortran.dg/bessel_5.f90 new file mode 100644 index 000000000000..aab45cafe17f --- /dev/null +++ b/gcc/testsuite/gfortran.dg/bessel_5.f90 @@ -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