Accurate numerical evaluation of a ratio of Bessel functions For fixed positive $x$, I would like to accurately (close to full precision) evaluate the ratio of Bessel functions of the first kind
$$
R_n(x):= \frac{J_{n+1}(x)}{J_n(x)}
$$
as $n$ becomes extremely large. Evaluating the special functions directly fails because both numerator and denominator become very small. However, it is known that this ratio obeys the large $n$ asymptotics
$$
R_n(x) \sim \sqrt{\frac{n}{n+1}} \frac{x}{2(n+1)}
$$
so in principle it seems possible that there is a stable method for computing this quantity. Does anybody have any suggestions?
 A: Using 
$$J_{\nu}(x) \approx \frac{1}{\sqrt{2 \, \pi \, \nu}} \, \left(\frac{e \, x}{2 \, \nu}\right)^{n}$$
for $\nu \to \infty$ leads to
$$R_{n}(x) = \frac{J_{n+1}(x)}{J_{n}(x)} \approx \frac{e \, x}{2} \, \frac{1}{n+1} \, \left(\frac{n}{n+1}\right)^{n+ \frac{1}{2}}.$$
If this is not within the realm sought then continued fractions, approximation of integral representations, etc may be applied. Similar results may be obtained though.
A: Boost.math uses Miller's algorithm for this problem; see J.M.'s description here. Code to compute the ratio to 100 digits or so:
#include <iostream>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

using std::sqrt;
using boost::math::cyl_bessel_j;
using namespace boost::multiprecision;

int main(int argc, char** argv)
{
    cpp_dec_float_100 x{5.7};
    for (int i = 10000; i < 50000; ++i)
    {
        std::cout << "J_" << i << "/J_" << i+1 << " = " << cyl_bessel_j<cpp_dec_float_100>(i+1, x)/cyl_bessel_j<cpp_dec_float_100>(i, x) << "\n";
        std::cout << "asymptotic form = " << sqrt((cpp_dec_float_100)i /(cpp_dec_float_100) i+1)*x/(2*i+2) << "\n";
    }
}

