Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders.
Summary: cyl_bessel_j returns wrong result for x>1000 for high orders.
Status: RESOLVED FIXED
Alias: None
Product: gcc
Classification: Unclassified
Component: libstdc++ (show other bugs)
Version: 7.2.0
: P3 normal
Target Milestone: ---
Assignee: Ed Smith-Rowland
URL:
Keywords: wrong-code
Depends on:
Blocks:
 
Reported: 2017-12-23 23:02 UTC by Ruslan
Modified: 2018-11-21 21:57 UTC (History)
3 users (show)

See Also:
Host:
Target:
Build:
Known to work:
Known to fail:
Last reconfirmed: 2018-01-04 00:00:00

Attachments
Add an attachment (proposed patch, testcase, etc.)

Note You need to log in before you can comment on or make changes to this bug.
Description Ruslan 2017-12-23 23:02:05 UTC 
The following test program prints results of C++17 std::cyl_bessel_j(100,1000.0001) and corresponding result given by GSL:

#include <cmath>
#include <iostream>
#include <gsl/gsl_sf_bessel.h>

int main()
{
    const double volatile n = 100;
    const double volatile x = 1000.0001;
    std::cout.precision(std::numeric_limits<decltype(x)>::digits10);
    const auto valueCXX17 = std::cyl_bessel_j(n,x);
    const auto valueGSL   = gsl_sf_bessel_Jn (n,x);
    std::cout << "C++17: " << valueCXX17 << "\n"
              << "GSL  : " << valueGSL << "\n";
}

I get the following output:

C++17: 0.433818396252946
GSL  : 0.0116783669817645

Comparison with Boost.Math and Wolfram Mathematica shows that GSL is right, while stdc++ is wrong.

For x<=1000 there's no such problem. As n decreases, the imprecision gradually gets smaller.
Comment 1 Ruslan 2017-12-23 23:04:28 UTC 
> As n decreases, the imprecision gradually gets smaller.
To avoid confusion: this statement is for fixed x>1000.
Comment 2 Michele Pezzutti 2018-01-01 02:39:58 UTC 
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.
Patch proposed to mailing list.
Comment 3 Ed Smith-Rowland 2018-11-18 18:32:58 UTC 
Author: emsr
Date: Sun Nov 18 18:32:26 2018
New Revision: 266252

URL: https://gcc.gnu.org/viewcvs?rev=266252&root=gcc&view=rev
Log:
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.


Modified:
    trunk/libstdc++-v3/ChangeLog
    trunk/libstdc++-v3/include/tr1/bessel_function.tcc
    trunk/libstdc++-v3/testsuite/special_functions/08_cyl_bessel_j/check_value.cc
    trunk/libstdc++-v3/testsuite/special_functions/10_cyl_neumann/check_value.cc
    trunk/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/09_cyl_bessel_j/check_value.cc
    trunk/libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/11_cyl_neumann/check_value.cc
Comment 4 Ed Smith-Rowland 2018-11-21 21:57:21 UTC 
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.