Re: [Patch, fortran] PR83076 - [8 Regression] ICE in gfc_deallocate_scalar_with_status, at fortran/trans.c:1598
Committed to trunk as revision 256065. Damian, it would be good if you would confirm that there are no issues with applying the patch to 7-branch. Thanks for all the help. Paul On 28 December 2017 at 19:58, Damian Rouson wrote: > I applied the patch the trunk and confirmed that it doesn’t break any > previously > passing OpenCoarrays tests. Is that sufficient or should I also try applying > the > patch to the 7 branch? > > Damian > > -- "If you can't explain it simply, you don't understand it well enough" - Albert Einstein
[PATCH] avoid -Wsuggest-attribute=const for void functions (PR 83559)
Among the recent -Wattributes enhancements is to diagnose declarations of void functions with attribute const or pure. Declaring a void function const or pure means that calls to the function have no effect and could be (and, with optimization, are) eliminated. Thus the attribute on such a function would imply that either the function is useless or the attribute is (when the function is called inidirectly). As it happens, however, -Wsuggest-attribute=const suggests adding the const attribute to candidate functions that return void, thus causing a conflict between the two options. The attached patch resolves the conflict by avoiding -Wsuggest-attribute=const for void functions. Martin PR c/83559 - -Wsuggest-attribute=const conflicts with -Wattributes warning about const attribute on function returning void gcc/ChangeLog: PR c/83559 * doc/extend.texi (attribute const): Fix a typo. * ipa-pure-const.c ((warn_function_const, warn_function_pure): Avoid issuing -Wsuggest-attribute for void functions. gcc/testsuite/ChangeLog: PR c/83559 * gcc.dg/const-2.c: New test. * gcc.dg/pure-3.c: New test. diff --git a/gcc/doc/extend.texi b/gcc/doc/extend.texi index 2a553ad..b2339f6 100644 --- a/gcc/doc/extend.texi +++ b/gcc/doc/extend.texi @@ -2500,7 +2500,7 @@ definition than the similar @code{pure} attribute below because it prohibits the function from reading global variables. Consequently, the presence of the attribute on a function declarations allows GCC to emit more efficient code for some calls to the function. Decorating the same function with -both the @code{const} and the @code{pure} attribute is diagnnosed. +both the @code{const} and the @code{pure} attribute is diagnosed. @cindex pointer arguments Note that a function that has pointer arguments and examines the data diff --git a/gcc/ipa-pure-const.c b/gcc/ipa-pure-const.c index 09ca359..f5eed37 100644 --- a/gcc/ipa-pure-const.c +++ b/gcc/ipa-pure-const.c @@ -213,9 +213,13 @@ suggest_attribute (int option, tree decl, bool known_finite, static void warn_function_pure (tree decl, bool known_finite) { - static hash_set *warned_about; + /* Declaring a void function pure makes no sense and is diagnosed + by -Wattributes because calling it would have no effect. */ + if (VOID_TYPE_P (TREE_TYPE (TREE_TYPE (decl +return; - warned_about + static hash_set *warned_about; + warned_about = suggest_attribute (OPT_Wsuggest_attribute_pure, decl, known_finite, warned_about, "pure"); } @@ -226,8 +230,13 @@ warn_function_pure (tree decl, bool known_finite) static void warn_function_const (tree decl, bool known_finite) { + /* Declaring a void function const makes no sense is diagnosed + by -Wattributes because calling it would have no effect. */ + if (VOID_TYPE_P (TREE_TYPE (TREE_TYPE (decl +return; + static hash_set *warned_about; - warned_about + warned_about = suggest_attribute (OPT_Wsuggest_attribute_const, decl, known_finite, warned_about, "const"); } diff --git a/gcc/testsuite/gcc.dg/const-2.c b/gcc/testsuite/gcc.dg/const-2.c new file mode 100644 index 000..e48005d --- /dev/null +++ b/gcc/testsuite/gcc.dg/const-2.c @@ -0,0 +1,22 @@ +/* PR tree-optimization/83559 - -Wsuggest-attribute=const conflicts with + -Wattributes warning about const attribute on function returning void + { dg-do compile { target nonpic } } + { dg-options "-O2 -Wsuggest-attribute=const" } */ + +int f_i_v (void)/* { dg-warning "candidate for attribute .const." } */ +{ + return 0; +} + +int f_i () /* { dg-warning "candidate for attribute .const." } */ +{ + return 0; +} + +void f_v_v (void) /* { dg-bogus "candidate" } */ +{ +} + +void f_v () /* { dg-bogus "candidate" } */ +{ +} diff --git a/gcc/testsuite/gcc.dg/pure-3.c b/gcc/testsuite/gcc.dg/pure-3.c new file mode 100644 index 000..2eeb8a3 --- /dev/null +++ b/gcc/testsuite/gcc.dg/pure-3.c @@ -0,0 +1,24 @@ +/* PR tree-optimization/83559 - -Wsuggest-attribute=const conflicts with + -Wattributes warning about const attribute on function returning void + { dg-do compile { target nonpic } } + { dg-options "-O2 -Wsuggest-attribute=pure" } */ + +int global; + +int f_i_v (void)/* { dg-warning "candidate for attribute .pure." } */ +{ + return global; +} + +int f_i () /* { dg-warning "candidate for attribute .pure." } */ +{ + return global; +} + +void f_v_v (void) /* { dg-bogus "candidate" } */ +{ +} + +void f_v () /* { dg-bogus "candidate" } */ +{ +}
[PATCH] handle invalid calls to built-ins with no prototype (PR 83603)
The -Wrestrict code assumes that built-ins are called with the correct number of arguments. When this isn't so it crashes. The attached patch avoids the ICE due to this error. There are outstanding assumptions that the type of actual arguments to built-ins declared without a prototype matches the expected type. Those will be corrected in a subsequent patch. Martin PR tree-optimization/83603 - ICE in builtin_memref at gcc/gimple-ssa-warn-restrict.c:238 gcc/ChangeLog: PR tree-optimization/83603 * calls.c (maybe_warn_nonstring_arg): Avoid accessing function arguments past the endof the argument list in functions declared without a prototype. * gimple-ssa-warn-restrict.c (wrestrict_dom_walker::check_call): Avoid checking when arguments are null. diff --git a/gcc/calls.c b/gcc/calls.c index 9b7e118..0cfc20a 100644 --- a/gcc/calls.c +++ b/gcc/calls.c @@ -1606,6 +1606,8 @@ maybe_warn_nonstring_arg (tree fndecl, tree exp) bool with_bounds = CALL_WITH_BOUNDS_P (exp); + unsigned nargs = call_expr_nargs (exp); + /* The bound argument to a bounded string function like strncpy. */ tree bound = NULL_TREE; @@ -1620,12 +1622,20 @@ maybe_warn_nonstring_arg (tree fndecl, tree exp) case BUILT_IN_STRNCASECMP: case BUILT_IN_STRNCPY: case BUILT_IN_STRNCPY_CHK: - bound = CALL_EXPR_ARG (exp, with_bounds ? 4 : 2); - break; + { + unsigned argno = with_bounds ? 4 : 2; + if (argno < nargs) + bound = CALL_EXPR_ARG (exp, argno); + break; + } case BUILT_IN_STRNDUP: - bound = CALL_EXPR_ARG (exp, with_bounds ? 2 : 1); - break; + { + unsigned argno = with_bounds ? 2 : 1; + if (argno < nargs) + bound = CALL_EXPR_ARG (exp, argno); + break; + } default: break; @@ -1645,6 +1655,11 @@ maybe_warn_nonstring_arg (tree fndecl, tree exp) for (unsigned argno = 0; ; ++argno, function_args_iter_next (&it)) { + /* Avoid iterating past the declared argument in a call + to function declared without a prototype. */ + if (argno >= nargs) + break; + tree argtype = function_args_iter_cond (&it); if (!argtype) break; diff --git a/gcc/gimple-ssa-warn-restrict.c b/gcc/gimple-ssa-warn-restrict.c index ebec381..a2a29cd 100644 --- a/gcc/gimple-ssa-warn-restrict.c +++ b/gcc/gimple-ssa-warn-restrict.c @@ -1710,7 +1710,9 @@ wrestrict_dom_walker::check_call (gcall *call) if (!dstwr && strfun) dstwr = size_one_node; - if (check_bounds_or_overlap (call, dst, src, dstwr, NULL_TREE)) + /* DST and SRC can be null for a call with an insufficient number + of arguments to a built-in function declared without a protype. */ + if (!dst || !src || check_bounds_or_overlap (call, dst, src, dstwr, NULL_TREE)) return; /* Avoid diagnosing the call again. */ diff --git a/gcc/testsuite/gcc.dg/Wrestrict-4.c b/gcc/testsuite/gcc.dg/Wrestrict-4.c new file mode 100644 index 000..f2398ef --- /dev/null +++ b/gcc/testsuite/gcc.dg/Wrestrict-4.c @@ -0,0 +1,110 @@ +/* PR tree-optimization/83603 - ICE in builtin_memref at + gcc/gimple-ssa-warn-restrict.c:238 + Test to verify that invalid calls to built-in functions declared + without a prototype don't cause an ICE. + { dg-do compile } + { dg-options "-O2 -Warray-bounds -Wrestrict" } */ + +void* memcpy (); +void* memmove (); +char* stpcpy (); +char* strcat (); +char* strcpy (); +char* strncat (); +char* strncpy (); + +void* test_memcpy_0 () +{ + return memcpy (); +} + +void* test_memcpy_1 (void *d) +{ + return memcpy (d); +} + +void* test_memcpy_2 (void *d, const void *s) +{ + return memcpy (d, s); +} + + +void* test_memmove_0 () +{ + return memmove (); +} + +void* test_memmove_1 (void *d) +{ + return memmove (d); +} + +void* test_memmove_2 (void *d, const void *s) +{ + return memmove (d, s); +} + + +void* test_stpcpy_0 () +{ + return stpcpy (); +} + +void* test_stpcpy_1 (char *d) +{ + return stpcpy (d); +} + + +char* test_strcat_0 () +{ + return strcat (); +} + +char* test_strcat_1 (char *d) +{ + return strcat (d); +} + + +void* test_strcpy_0 () +{ + return strcpy (); +} + +void* test_strcpy_1 (char *d) +{ + return strcpy (d); +} + + +char* test_strncat_0 () +{ + return strncat (); +} + +char* test_strncat_1 (char *d) +{ + return strncat (d); +} + +char* test_strncat_2 (char *d, const char *s) +{ + return strncat (d, s); +} + + +void* test_strncpy_0 () +{ + return strncpy (); +} + +void* test_strncpy_1 (char *d) +{ + return strncpy (d); +} + +void* test_strncpy_2 (char *d, const char *s) +{ + return strncpy (d, s); +}
[PATCH] constrain strcat destination offset (PR 83642)
The ICE in the test case submitted in PR tree-optimization/83640 is triggered by the tree-ssa-strlen pass transforming calls to strcat to strcpy with an offset pointing to the terminating NUL of the destination string, and allowing the upper bound of the offset's range to exceed PTRDIFF_MAX by the length of the appended constant string. Although the ICE itself is due to an unsafe assumption in the -Wrestrict checker, the excessive upper bound in the strcat case suggests an optimization opportunity that's missing from the tree-ssa-strlen pass: namely, to reduce the offset's upper bound by the length of the appended string. The attached patch adds this minor optimization. Martin PR tree-optimization/83642 - excessive strlen range after a strcat of non-empty string gcc/ChangeLog: PR tree-optimization/83642 * tree-ssa-strlen.c (handle_builtin_strcat): Set offset range. gcc/testsuite/ChangeLog: PR tree-optimization/83642 * gcc.dg/strlenopt-39.c: New test. diff --git a/gcc/testsuite/gcc.dg/strlenopt-39.c b/gcc/testsuite/gcc.dg/strlenopt-39.c new file mode 100644 index 000..b070b6b --- /dev/null +++ b/gcc/testsuite/gcc.dg/strlenopt-39.c @@ -0,0 +1,84 @@ +/* PR tree-optimization/83642] excessive strlen range after a strcat + of non-empty string + { dg-do compile } + { dg-options "-O2 -fdump-tree-optimized" } */ + +void test_strcat_0_strlen_range (char *dst, char *src) +{ + __SIZE_TYPE__ n = __builtin_strlen (dst); + + __builtin_strcat (dst, ""); + __builtin_strcat (dst, src); + + if (n >= __PTRDIFF_MAX__) +__builtin_abort (); +} + +void test_strcat_1_strlen_range (char *dst, char *src) +{ + __SIZE_TYPE__ n = __builtin_strlen (dst); + + __builtin_strcat (dst, "1"); + __builtin_strcat (dst, src); + + if (n >= __PTRDIFF_MAX__ - 1) +__builtin_abort (); +} + +void test_strcat_2_strlen_range (char *dst, char *src) +{ + __SIZE_TYPE__ n = __builtin_strlen (dst); + + __builtin_strcat (dst, "12"); + __builtin_strcat (dst, src); + + if (n >= __PTRDIFF_MAX__ - 2) +__builtin_abort (); +} + +void test_strcat_3_strlen_range (char *dst, char *src) +{ + __SIZE_TYPE__ n = __builtin_strlen (dst); + char a[] = "123"; + + __builtin_strcat (dst, a); + __builtin_strcat (dst, src); + + if (n >= __PTRDIFF_MAX__ - 3) +__builtin_abort (); +} + +void test_stpcpy_7_strlen_range (char *dst, char *src) +{ + __SIZE_TYPE__ n = __builtin_strlen (dst); + + __builtin_stpcpy (dst + n, "1234567"); + __builtin_strcat (dst, src); + + if (n >= __PTRDIFF_MAX__ - 7) +__builtin_abort (); +} + +void test_strcpy_8_strlen_range (char *dst, char *src) +{ + __SIZE_TYPE__ n = __builtin_strlen (dst); + + __builtin_strcpy (dst + n, "12345678"); + __builtin_strcat (dst, src); + + if (n >= __PTRDIFF_MAX__ - 8) +__builtin_abort (); +} + +void test_memcpy_9_strlen_range (char *dst, char *src) +{ + __SIZE_TYPE__ n = __builtin_strlen (dst); + + __builtin_memcpy (dst + n, "123456789", 10); + __builtin_strcat (dst, src); + + if (n >= __PTRDIFF_MAX__ - 9) +__builtin_abort (); +} + +/* { dg-final { scan-tree-dump-not "abort \\(\\)" "optimized" } } */ diff --git a/gcc/tree-ssa-strlen.c b/gcc/tree-ssa-strlen.c index be6ab9f..9e42232 100644 --- a/gcc/tree-ssa-strlen.c +++ b/gcc/tree-ssa-strlen.c @@ -2486,10 +2486,35 @@ handle_builtin_strcat (enum built_in_function bcode, gimple_stmt_iterator *gsi) if (endptr) dst = fold_convert_loc (loc, TREE_TYPE (dst), unshare_expr (endptr)); else -dst = fold_build2_loc (loc, POINTER_PLUS_EXPR, - TREE_TYPE (dst), unshare_expr (dst), - fold_convert_loc (loc, sizetype, - unshare_expr (dstlen))); +{ + if (TREE_CODE (dstlen) == PLUS_EXPR + && TREE_CODE (TREE_OPERAND (dstlen, 0)) == SSA_NAME + && TREE_CODE (TREE_OPERAND (dstlen, 1)) == INTEGER_CST) + { + /* For a call to strcat (DST, SRC) where DST = D + OFFSET and + where OFFSET = VAROFF + CSTOFF, adjust the upper bound of + VAROFF's range by the constant CSTOFF so that the result + is still guaranteed to be less than PTRDIFF_MAX, the size + of the longest possible string. */ + tree varoff = TREE_OPERAND (dstlen, 0); + wide_int minoff, maxoff; + value_range_type rng = get_range_info (varoff, &minoff, &maxoff); + + if (rng == VR_RANGE) + { + tree cstoff = TREE_OPERAND (dstlen, 1); + + maxoff -= wi::to_wide (cstoff); + set_range_info (varoff, rng, minoff, maxoff); + } + } + + dst = fold_build2_loc (loc, POINTER_PLUS_EXPR, + TREE_TYPE (dst), unshare_expr (dst), + fold_convert_loc (loc, sizetype, + unshare_expr (dstlen))); +} + dst = force_gimple_operand_gsi (gsi, dst, true, NULL_TREE, true, GSI_SAME_STMT); if (dump_file && (dump_flags & TDF_DETAILS) != 0)
[PATCH] correct handling of offset ranges that cross PTRDIFF_MAX (PR 83640)
PR tree-optimization/83640 - ice in generic_overlap at gimple-ssa-warn-restrict.c:814 highlights out a class of cases where the -Wrestrict checker assumes that the range of a pointer offset that us represented by a VR_RANGE has a lower bound that is less than its upper bound. The pass asserts that this is so and cases to the contrary trigger an ICE. The submitted test case that triggers this ICE is due to a missed optimization in tree-ssa-strlen (which could, and IMO should, guarantee that the offsets it constructs are in such a range, and I'll submit a separate patch with that change), but the assumption isn't safe in general for offsets whose range happens to straddle the PTRDIFF_MAX boundary, i.e., whose lower bound is greater than its upper bound. The attached patch makes adjustments to remove this assumption and avoid the ICE. Tested on x86_64-linux with no regressions. Martin PS The patch also removes the same troublesome assertion that's also removed in the one below: https://gcc.gnu.org/ml/gcc-patches/2017-12/msg01399.html PR tree-optimization/83640 - ice in generic_overlap at gimple-ssa-warn-restrict.c:814 gcc/ChangeLog: PR tree-optimization/83640 * gimple-ssa-warn-restrict.c (builtin_access::builtin_access): Avoid subtracting negative offset from size. (builtin_access::overlap): Adjust offset bounds of the access to fall within the size of the object if possible. (maybe_diag_overlap): Remove troublesome assertion. gcc/testsuite/ChangeLog: PR tree-optimization/83640 * gcc.dg/Wrestrict-3.c: New test. * gcc.dg/pr83640.c: New test. diff --git a/gcc/gimple-ssa-warn-restrict.c b/gcc/gimple-ssa-warn-restrict.c index ac545e4..ebec381 100644 --- a/gcc/gimple-ssa-warn-restrict.c +++ b/gcc/gimple-ssa-warn-restrict.c @@ -698,9 +698,10 @@ builtin_access::builtin_access (gcall *call, builtin_memref &dst, /* For string functions, adjust the size range of the source reference by the inverse boundaries of the offset (because - the higher the offset into the string the shorter its + the higher the offset into the string the shorter its length). */ - if (srcref->offrange[1] < srcref->sizrange[0]) + if (srcref->offrange[1] >= 0 + && srcref->offrange[1] < srcref->sizrange[0]) srcref->sizrange[0] -= srcref->offrange[1]; else srcref->sizrange[0] = 0; @@ -1134,30 +1135,53 @@ builtin_access::overlap () if (!dstref->base || !srcref->base) return false; - /* If the base object is an array adjust the lower bound of the offset - to be non-negative. */ + /* Set the access offsets. */ + acs.dstoff[0] = dstref->offrange[0]; + acs.dstoff[1] = dstref->offrange[1]; + + /* If the base object is an array adjust the bounds of the offset + to be non-negative and within the bounds of the array if possible. */ if (dstref->base && TREE_CODE (TREE_TYPE (dstref->base)) == ARRAY_TYPE) -acs.dstoff[0] = wi::smax (dstref->offrange[0], 0); - else -acs.dstoff[0] = dstref->offrange[0]; +{ + if (acs.dstoff[0] < 0 && acs.dstoff[1] >= 0) + acs.dstoff[0] = 0; - acs.dstoff[1] = dstref->offrange[1]; + if (acs.dstoff[1] < acs.dstoff[0]) + { + if (tree size = TYPE_SIZE_UNIT (TREE_TYPE (dstref->base))) + acs.dstoff[1] = wi::umin (acs.dstoff[1], wi::to_offset (size)); + else + acs.dstoff[1] = wi::umin (acs.dstoff[1], maxobjsize); + } +} + + acs.srcoff[0] = srcref->offrange[0]; + acs.srcoff[1] = srcref->offrange[1]; if (srcref->base && TREE_CODE (TREE_TYPE (srcref->base)) == ARRAY_TYPE) -acs.srcoff[0] = wi::smax (srcref->offrange[0], 0); - else -acs.srcoff[0] = srcref->offrange[0]; +{ + if (acs.srcoff[0] < 0 && acs.srcoff[1] >= 0) + acs.srcoff[0] = 0; - acs.srcoff[1] = srcref->offrange[1]; + if (tree size = TYPE_SIZE_UNIT (TREE_TYPE (srcref->base))) + acs.srcoff[1] = wi::umin (acs.srcoff[1], wi::to_offset (size)); + else if (acs.srcoff[1] < acs.srcoff[0]) + acs.srcoff[1] = wi::umin (acs.srcoff[1], maxobjsize); +} - /* When the lower bound of the offset is less that the upper bound - disregard it and use the inverse of the maximum object size - instead. The upper bound is the result of a negative offset - being represented as a large positive value. */ + /* When the upper bound of the offset is less than the lower bound + the former is the result of a negative offset being represented + as a large positive value or vice versa. The resulting range is + a union of two subranges: [MIN, UB] and [LB, MAX]. Since such + a union is not representable using the current data structure + replace it with the full range of offsets. */ if (acs.dstoff[1] < acs.dstoff[0]) -acs.dstoff[0] = -maxobjsize; +{ + acs.dstoff[0] = -maxobjsize - 1; + acs.dstoff[1] = maxobjsize; +} /* Validate the offset and size of each reference on its own first. This is independent of whether or not the base objects a
Re: Fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders
On 12/31/2017 09:38 PM, Michele Pezzutti wrote: Hi. This patch intends to fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders. Seehttps://gcc.gnu.org/bugzilla/show_bug.cgi?id=83566 forissue description. * libstdc++-v3/include/tr1/bessel_function.tcc Series expansion in __cyl_bessel_jn_asymp() shall not be truncated at the first terms. At least nu/2 terms shall be added, in order to have a guaranteed error bound. * libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc Testcase for x > 1000 added. * libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc Testcase for x > 1000 added. diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7ac733d..852a31e 100644 --- a/libstdc++-v3/include/tr1/bessel_function.tcc +++ b/libstdc++-v3/include/tr1/bessel_function.tcc @@ -359,15 +359,32 @@ namespace tr1 __cyl_bessel_jn_asymp(_Tp __nu, _Tp __x, _Tp & __Jnu, _Tp & __Nnu) { const _Tp __mu = _Tp(4) * __nu * __nu; - const _Tp __mum1 = __mu - _Tp(1); - const _Tp __mum9 = __mu - _Tp(9); - const _Tp __mum25 = __mu - _Tp(25); - const _Tp __mum49 = __mu - _Tp(49); - const _Tp __xx = _Tp(64) * __x * __x; - const _Tp __P = _Tp(1) - __mum1 * __mum9 / (_Tp(2) * __xx) - * (_Tp(1) - __mum25 * __mum49 / (_Tp(12) * __xx)); - const _Tp __Q = __mum1 / (_Tp(8) * __x) - * (_Tp(1) - __mum9 * __mum25 / (_Tp(6) * __xx)); + const _Tp __8x = _Tp(8) * __x; + + _Tp __P = _Tp(1); + _Tp __Q = _Tp(0); + + _Tp __eps = std::numeric_limits<_Tp>::epsilon(); + + _Tp __term = _Tp(1); + _Tp __epsP = _Tp(0); + _Tp __epsQ = _Tp(0); + + unsigned long __2k; + for (unsigned long k = 1; k < 1000; k+=2) { + __2k = 2 * k; + + __term *= (__mu - _Tp((__2k - 1) * (__2k - 1))) / (k * __8x); + __Q += __term; + __epsQ = std::abs(__term) < std::abs(__eps * __Q); + + __term *= -(__mu - _Tp((__2k + 1) * (__2k + 1)))/ (_Tp(k + 1) * __8x); + __P += __term; + __epsP = std::abs(__term) < std::abs(__eps * __P); + + if (__epsP && __epsQ && k > __nu / 2.) + break; + } const _Tp __chi = __x - (__nu + _Tp(0.5L)) * __numeric_constants<_Tp>::__pi_2(); diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc index 26f4dd3..e340b78 100644 --- a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc @@ -698,6 +698,26 @@ data026[21] = }; const double toler026 = 1.0006e-11; +// Test data for nu=100.00. +// max(|f - f_GSL|): 2.5857788132910287e-14 +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 +const testcase_cyl_bessel_j +data027[11] = +{ + { 0.0116761350077845, 100., 1000. }, + {-0.0116998547780258, 100., 1100. }, + {-0.0228014834050837, 100., 1200. }, + {-0.0169735007873739, 100., 1300. }, + {-0.0014154528803530, 100., 1400. }, + { 0.017265844988, 100., 1500. }, + { 0.0198025620201474, 100., 1600. }, + { 0.0161297712798388, 100., 1700. }, + { 0.0053753369281577, 100., 1800. }, + {-0.0069238868725646, 100., 1900. }, + {-0.0154878717200738, 100., 2000. }, +}; +const double toler027 = 1.0006e-10; + template void test(const testcase_cyl_bessel_j (&data)[Num], Ret toler) @@ -748,5 +768,6 @@ main() test(data024, toler024); test(data025, toler025); test(data026, toler026); + test(data027, toler027); return 0; } diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc index 5579149..9caf836 100644 --- a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc @@ -742,6 +742,26 @@ data028[20] = }; const double toler028 = 1.0006e-11; +// Test data for nu=100.00. +// max(|f - f_GSL|): 3.1049815496508870e-14 +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e
Re: Fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders
On 12/31/2017 09:38 PM, Michele Pezzutti wrote: Hi. This patch intends to fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders. Seehttps://gcc.gnu.org/bugzilla/show_bug.cgi?id=83566 forissue description. * libstdc++-v3/include/tr1/bessel_function.tcc Series expansion in __cyl_bessel_jn_asymp() shall not be truncated at the first terms. At least nu/2 terms shall be added, in order to have a guaranteed error bound. * libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc Testcase for x > 1000 added. * libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc Testcase for x > 1000 added. diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc index 7ac733d..852a31e 100644 --- a/libstdc++-v3/include/tr1/bessel_function.tcc +++ b/libstdc++-v3/include/tr1/bessel_function.tcc @@ -359,15 +359,32 @@ namespace tr1 __cyl_bessel_jn_asymp(_Tp __nu, _Tp __x, _Tp & __Jnu, _Tp & __Nnu) { const _Tp __mu = _Tp(4) * __nu * __nu; - const _Tp __mum1 = __mu - _Tp(1); - const _Tp __mum9 = __mu - _Tp(9); - const _Tp __mum25 = __mu - _Tp(25); - const _Tp __mum49 = __mu - _Tp(49); - const _Tp __xx = _Tp(64) * __x * __x; - const _Tp __P = _Tp(1) - __mum1 * __mum9 / (_Tp(2) * __xx) - * (_Tp(1) - __mum25 * __mum49 / (_Tp(12) * __xx)); - const _Tp __Q = __mum1 / (_Tp(8) * __x) - * (_Tp(1) - __mum9 * __mum25 / (_Tp(6) * __xx)); + const _Tp __8x = _Tp(8) * __x; + + _Tp __P = _Tp(1); + _Tp __Q = _Tp(0); + + _Tp __eps = std::numeric_limits<_Tp>::epsilon(); + + _Tp __term = _Tp(1); + _Tp __epsP = _Tp(0); + _Tp __epsQ = _Tp(0); + + unsigned long __2k; + for (unsigned long k = 1; k < 1000; k+=2) { + __2k = 2 * k; + + __term *= (__mu - _Tp((__2k - 1) * (__2k - 1))) / (k * __8x); + __Q += __term; + __epsQ = std::abs(__term) < std::abs(__eps * __Q); + + __term *= -(__mu - _Tp((__2k + 1) * (__2k + 1)))/ (_Tp(k + 1) * __8x); + __P += __term; + __epsP = std::abs(__term) < std::abs(__eps * __P); + + if (__epsP && __epsQ && k > __nu / 2.) + break; + } const _Tp __chi = __x - (__nu + _Tp(0.5L)) * __numeric_constants<_Tp>::__pi_2(); diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc index 26f4dd3..e340b78 100644 --- a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc @@ -698,6 +698,26 @@ data026[21] = }; const double toler026 = 1.0006e-11; +// Test data for nu=100.00. +// max(|f - f_GSL|): 2.5857788132910287e-14 +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e-11 +const testcase_cyl_bessel_j +data027[11] = +{ + { 0.0116761350077845, 100., 1000. }, + {-0.0116998547780258, 100., 1100. }, + {-0.0228014834050837, 100., 1200. }, + {-0.0169735007873739, 100., 1300. }, + {-0.0014154528803530, 100., 1400. }, + { 0.017265844988, 100., 1500. }, + { 0.0198025620201474, 100., 1600. }, + { 0.0161297712798388, 100., 1700. }, + { 0.0053753369281577, 100., 1800. }, + {-0.0069238868725646, 100., 1900. }, + {-0.0154878717200738, 100., 2000. }, +}; +const double toler027 = 1.0006e-10; + template void test(const testcase_cyl_bessel_j (&data)[Num], Ret toler) @@ -748,5 +768,6 @@ main() test(data024, toler024); test(data025, toler025); test(data026, toler026); + test(data027, toler027); return 0; } diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc index 5579149..9caf836 100644 --- a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc @@ -742,6 +742,26 @@ data028[20] = }; const double toler028 = 1.0006e-11; +// Test data for nu=100.00. +// max(|f - f_GSL|): 3.1049815496508870e-14 +// max(|f - f_GSL| / |f_GSL|): 1.6767662425535933e