PR libstdc++/83566 - cyl_bessel_j returns wrong result for x>1000

2018-11-16  Michele Pezzutti <mpezz@tiscali.it>
	    Edward Smith-Rowland  <3dw4rd@verizon.net>

	PR libstdc++/83566 - cyl_bessel_j returns wrong result for x>1000
	for high orders.
	* include/tr1/bessel_function.tcc: Perform no fewer than nu/2 iterations
	of the asymptotic series (nu is the Bessel order).
	* testsuite/tr1/5_numerical_facilities/special_functions/
	09_cyl_bessel_j/check_value.cc: Add tests at nu=100, 1000<=x<=2000.
	* testsuite/tr1/5_numerical_facilities/special_functions/	
	11_cyl_neumann/check_value.cc: Ditto.
	* testsuite/special_functions/08_cyl_bessel_j/check_value.cc: Ditto.
	* testsuite/special_functions/10_cyl_neumann/check_value.cc: Ditto.


Co-Authored-By: Edward Smith-Rowland <3dw4rd@verizon.net>

From-SVN: r266252

2 files changed, 68 insertions, 0 deletions

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 06a2d19..ee587cb 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,39 @@ data026[21] =
 };
 const double toler026 = 1.0000000000000006e-11;
 
+// Test data for nu=100.0000000000000000
+// max(|f - f_GSL|): 3.9438938226332709e-14 at index 19
+// max(|f - f_GSL| / |f_GSL|): 2.0193411077170867e-11
+// mean(f - f_GSL): 1.6682360684660055e-15
+// variance(f - f_GSL): 5.3274331668346898e-28
+// stddev(f - f_GSL): 2.3081232997469372e-14
+const testcase_cyl_bessel_j<double>
+data027[21] =
+{
+  {  1.1676135007789573e-02, 100.0000000000000000, 1000.0000000000000000, 0.0 },
+  { -1.1699854778025796e-02, 100.0000000000000000, 1100.0000000000000000, 0.0 },
+  { -2.2801483405083697e-02, 100.0000000000000000, 1200.0000000000000000, 0.0 },
+  { -1.6973500787373915e-02, 100.0000000000000000, 1300.0000000000000000, 0.0 },
+  { -1.4154528803481308e-03, 100.0000000000000000, 1400.0000000000000000, 0.0 },
+  {  1.3333726584495232e-02, 100.0000000000000000, 1500.0000000000000000, 0.0 },
+  {  1.9802562020148559e-02, 100.0000000000000000, 1600.0000000000000000, 0.0 },
+  {  1.6129771279838816e-02, 100.0000000000000000, 1700.0000000000000000, 0.0 },
+  {  5.3753369281536031e-03, 100.0000000000000000, 1800.0000000000000000, 0.0 },
+  { -6.9238868725645785e-03, 100.0000000000000000, 1900.0000000000000000, 0.0 },
+  { -1.5487871720069789e-02, 100.0000000000000000, 2000.0000000000000000, 0.0 },
+  { -1.7275186717671070e-02, 100.0000000000000000, 2100.0000000000000000, 0.0 },
+  { -1.2233030525173150e-02, 100.0000000000000000, 2200.0000000000000000, 0.0 },
+  { -2.8518508672241900e-03, 100.0000000000000000, 2300.0000000000000000, 0.0 },
+  {  7.0784372270289329e-03, 100.0000000000000000, 2400.0000000000000000, 0.0 },
+  {  1.3955367586928166e-02, 100.0000000000000000, 2500.0000000000000000, 0.0 },
+  {  1.5574059842493392e-02, 100.0000000000000000, 2600.0000000000000000, 0.0 },
+  {  1.1718043044647556e-02, 100.0000000000000000, 2700.0000000000000000, 0.0 },
+  {  4.0320953231285607e-03, 100.0000000000000000, 2800.0000000000000000, 0.0 },
+  { -4.6895111783053977e-03, 100.0000000000000000, 2900.0000000000000000, 0.0 },
+  { -1.1507715400035966e-02, 100.0000000000000000, 3000.0000000000000000, 0.0 },
+};
+const double toler027 = 1.0000000000000006e-10;
+
 template<typename Ret, unsigned int Num>
   void
   test(const testcase_cyl_bessel_j<Ret> (&data)[Num], Ret toler)
@@ -748,5 +781,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 3c0cbe7..d73f6a6 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,39 @@ data028[20] =
 };
 const double toler028 = 1.0000000000000006e-11;
 
+// Test data for nu=100.0000000000000000
+// max(|f - f_GSL|): 3.9022387751663778e-14 at index 16
+// max(|f - f_GSL| / |f_GSL|): 2.4760677072012703e-11
+// mean(f - f_GSL): 3.6878362466971231e-16
+// variance(f - f_GSL): 5.0707962306468580e-28
+// stddev(f - f_GSL): 2.2518428521206487e-14
+const testcase_cyl_neumann<double>
+data029[21] =
+{
+  { -2.2438688257729954e-02, 100.0000000000000000, 1000.0000000000000000, 0.0 },
+  { -2.1077595159819992e-02, 100.0000000000000000, 1100.0000000000000000, 0.0 },
+  { -3.5299439206692585e-03, 100.0000000000000000, 1200.0000000000000000, 0.0 },
+  {  1.4250019326536615e-02, 100.0000000000000000, 1300.0000000000000000, 0.0 },
+  {  2.1304679089735663e-02, 100.0000000000000000, 1400.0000000000000000, 0.0 },
+  {  1.5734395077905267e-02, 100.0000000000000000, 1500.0000000000000000, 0.0 },
+  {  2.5544633636137774e-03, 100.0000000000000000, 1600.0000000000000000, 0.0 },
+  { -1.0722045524849367e-02, 100.0000000000000000, 1700.0000000000000000, 0.0 },
+  { -1.8036919243226864e-02, 100.0000000000000000, 1800.0000000000000000, 0.0 },
+  { -1.6958415593079763e-02, 100.0000000000000000, 1900.0000000000000000, 0.0 },
+  { -8.8788704566276667e-03, 100.0000000000000000, 2000.0000000000000000, 0.0 },
+  {  2.2504407108413179e-03, 100.0000000000000000, 2100.0000000000000000, 0.0 },
+  {  1.1833215246712251e-02, 100.0000000000000000, 2200.0000000000000000, 0.0 },
+  {  1.6398784536343945e-02, 100.0000000000000000, 2300.0000000000000000, 0.0 },
+  {  1.4675984403642338e-02, 100.0000000000000000, 2400.0000000000000000, 0.0 },
+  {  7.7523920451654229e-03, 100.0000000000000000, 2500.0000000000000000, 0.0 },
+  { -1.5759822576003489e-03, 100.0000000000000000, 2600.0000000000000000, 0.0 },
+  { -9.9314877404787089e-03, 100.0000000000000000, 2700.0000000000000000, 0.0 },
+  { -1.4534495161704743e-02, 100.0000000000000000, 2800.0000000000000000, 0.0 },
+  { -1.4059273497237509e-02, 100.0000000000000000, 2900.0000000000000000, 0.0 },
+  { -8.9385158149605185e-03, 100.0000000000000000, 3000.0000000000000000, 0.0 },
+};
+const double toler029 = 1.0000000000000006e-10;
+
 template<typename Ret, unsigned int Num>
   void
   test(const testcase_cyl_neumann<Ret> (&data)[Num], Ret toler)
@@ -794,5 +827,6 @@ main()
   test(data026, toler026);
   test(data027, toler027);
   test(data028, toler028);
+  test(data029, toler029);
   return 0;
 }