Re: [Patch, fortran] PR83076 - [8 Regression] ICE in gfc_deallocate_scalar_with_status, at fortran/trans.c:1598

2018-01-01 Thread Paul Richard Thomas
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)

2018-01-01 Thread Martin Sebor

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)

2018-01-01 Thread Martin Sebor

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)

2018-01-01 Thread Martin Sebor

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)

2018-01-01 Thread Martin Sebor

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

2018-01-01 Thread Ed Smith-Rowland

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

2018-01-01 Thread Ed Smith-Rowland

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