From patchwork Tue Aug 17 21:14:23 2010 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Tobias Burnus X-Patchwork-Id: 61971 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 8169FB6F14 for ; Wed, 18 Aug 2010 07:14:45 +1000 (EST) Received: (qmail 9427 invoked by alias); 17 Aug 2010 21:14:40 -0000 Received: (qmail 9410 invoked by uid 22791); 17 Aug 2010 21:14:38 -0000 X-SWARE-Spam-Status: No, hits=-0.5 required=5.0 tests=AWL, BAYES_50, RCVD_IN_DNSWL_NONE X-Spam-Check-By: sourceware.org Received: from mx02.qsc.de (HELO mx02.qsc.de) (213.148.130.14) by sourceware.org (qpsmtpd/0.43rc1) with ESMTP; Tue, 17 Aug 2010 21:14:28 +0000 Received: from [192.168.178.22] (port-92-204-42-5.dynamic.qsc.de [92.204.42.5]) by mx02.qsc.de (Postfix) with ESMTP id 2310129600; Tue, 17 Aug 2010 23:14:23 +0200 (CEST) Message-ID: <4C6AFBAF.2030009@net-b.de> Date: Tue, 17 Aug 2010 23:14:23 +0200 From: Tobias Burnus User-Agent: Mozilla/5.0 (X11; U; Linux x86_64; de; 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> <4C6ABFE0.7030908@net-b.de> <20100817190154.GA36792@troutmask.apl.washington.edu> In-Reply-To: <20100817190154.GA36792@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 Steve Kargl wrote: > I'll look over the patch in more detail in a few hours. Thanks! Here is an updated version. With the help of Daniel K, I have found a solution to the constructor issue; JN now works as it should. YN also works, though for some of the examples, the result differs quite a bit from the result of mpfr_yn. I do not know whether I have chosen particularly difficult values or whether the algorithm has problems. Please check. (Cf. bessel_5.f90 - the ULP values given are the differences which maximally occur [though 1 ULP could be 0 ULP in some cases].) Build and currently regtesting on x86-64-linux. OK for the trunk? Tobias 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/intrinsic.c =================================================================== --- gcc/fortran/intrinsic.c (Revision 163318) +++ gcc/fortran/intrinsic.c (Arbeitskopie) @@ -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 163318) +++ gcc/fortran/intrinsic.h (Arbeitskopie) @@ -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/gfortran.h =================================================================== --- gcc/fortran/gfortran.h (Revision 163318) +++ gcc/fortran/gfortran.h (Arbeitskopie) @@ -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/check.c =================================================================== --- gcc/fortran/check.c (Revision 163318) +++ gcc/fortran/check.c (Arbeitskopie) @@ -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: gcc/fortran/intrinsic.texi =================================================================== --- gcc/fortran/intrinsic.texi (Revision 163318) +++ gcc/fortran/intrinsic.texi (Arbeitskopie) @@ -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 163318) +++ gcc/fortran/simplify.c (Arbeitskopie) @@ -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.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). */ + + 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.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, + -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 +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) { 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,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 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,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 ULP + +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 ULP + +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 ULP + +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 ULP + +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: 0 ULP + +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: 0 ULP + +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 ULP + +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 ULP + +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