[Bug libstdc++/83566] New: cyl_bessel_j returns wrong result for x>1000 for high orders.

[b7.10110111 at gmail dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=83566

            Bug ID: 83566
           Summary: cyl_bessel_j returns wrong result for x>1000 for high
                    orders.
           Product: gcc
           Version: 7.2.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libstdc++
          Assignee: unassigned at gcc dot gnu.org
          Reporter: b7.10110111 at gmail dot com
  Target Milestone: ---

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.