special_functions/detail/bessel_jy_asym.hpp

0001 //  Copyright (c) 2007 John Maddock
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 //
0007 // This is a partial header, do not include on it's own!!!
0008 //
0009 // Contains asymptotic expansions for Bessel J(v,x) and Y(v,x)
0010 // functions, as x -> INF.
0011 //
0012 #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
0013 #define BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
0014
0015 #ifdef _MSC_VER
0016 #pragma once
0017 #endif
0018
0019 #include <boost/math/special_functions/factorials.hpp>
0020
0021 namespace boost{ namespace math{ namespace detail{
0022
0023 template <class T>
0024 inline T asymptotic_bessel_amplitude(T v, T x)
0025 {
0026    // Calculate the amplitude of J(v, x) and Y(v, x) for large
0027    // x: see A&S 9.2.28.
0028    BOOST_MATH_STD_USING
0029    T s = 1;
0030    T mu = 4 * v * v;
0031    T txq = 2 * x;
0032    txq *= txq;
0033
0034    s += (mu - 1) / (2 * txq);
0035    s += 3 * (mu - 1) * (mu - 9) / (txq * txq * 8);
0036    s += 15 * (mu - 1) * (mu - 9) * (mu - 25) / (txq * txq * txq * 8 * 6);
0037
0038    return sqrt(s * 2 / (constants::pi<T>() * x));
0039 }
0040
0041 template <class T>
0042 T asymptotic_bessel_phase_mx(T v, T x)
0043 {
0044    //
0045    // Calculate the phase of J(v, x) and Y(v, x) for large x.
0046    // See A&S 9.2.29.
0047    // Note that the result returned is the phase less (x - PI(v/2 + 1/4))
0048    // which we'll factor in later when we calculate the sines/cosines of the result:
0049    //
0050    T mu = 4 * v * v;
0051    T denom = 4 * x;
0052    T denom_mult = denom * denom;
0053
0054    T s = 0;
0055    s += (mu - 1) / (2 * denom);
0056    denom *= denom_mult;
0057    s += (mu - 1) * (mu - 25) / (6 * denom);
0058    denom *= denom_mult;
0059    s += (mu - 1) * (mu * mu - 114 * mu + 1073) / (5 * denom);
0060    denom *= denom_mult;
0061    s += (mu - 1) * (5 * mu * mu * mu - 1535 * mu * mu + 54703 * mu - 375733) / (14 * denom);
0062    return s;
0063 }
0064
0065 template <class T, class Policy>
0066 inline T asymptotic_bessel_y_large_x_2(T v, T x, const Policy& pol)
0067 {
0068    // See A&S 9.2.19.
0069    BOOST_MATH_STD_USING
0070    // Get the phase and amplitude:
0071    T ampl = asymptotic_bessel_amplitude(v, x);
0072    T phase = asymptotic_bessel_phase_mx(v, x);
0073    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
0074    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
0075    //
0076    // Calculate the sine of the phase, using
0077    // sine/cosine addition rules to factor in
0078    // the x - PI(v/2 + 1/4) term not added to the
0079    // phase when we calculated it.
0080    //
0081    T cx = cos(x);
0082    T sx = sin(x);
0083    T ci = boost::math::cos_pi(v / 2 + 0.25f, pol);
0084    T si = boost::math::sin_pi(v / 2 + 0.25f, pol);
0085    T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
0086    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
0087    BOOST_MATH_INSTRUMENT_CODE(cos(x));
0088    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
0089    BOOST_MATH_INSTRUMENT_CODE(sin(x));
0090    return sin_phase * ampl;
0091 }
0092
0093 template <class T, class Policy>
0094 inline T asymptotic_bessel_j_large_x_2(T v, T x, const Policy& pol)
0095 {
0096    // See A&S 9.2.19.
0097    BOOST_MATH_STD_USING
0098    // Get the phase and amplitude:
0099    T ampl = asymptotic_bessel_amplitude(v, x);
0100    T phase = asymptotic_bessel_phase_mx(v, x);
0101    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
0102    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
0103    //
0104    // Calculate the sine of the phase, using
0105    // sine/cosine addition rules to factor in
0106    // the x - PI(v/2 + 1/4) term not added to the
0107    // phase when we calculated it.
0108    //
0109    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
0110    BOOST_MATH_INSTRUMENT_CODE(cos(x));
0111    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
0112    BOOST_MATH_INSTRUMENT_CODE(sin(x));
0113    T cx = cos(x);
0114    T sx = sin(x);
0115    T ci = boost::math::cos_pi(v / 2 + 0.25f, pol);
0116    T si = boost::math::sin_pi(v / 2 + 0.25f, pol);
0117    T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
0118    BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
0119    return sin_phase * ampl;
0120 }
0121
0122 template <class T>
0123 inline bool asymptotic_bessel_large_x_limit(int v, const T& x)
0124 {
0125    BOOST_MATH_STD_USING
0126       //
0127       // Determines if x is large enough compared to v to take the asymptotic
0128       // forms above.  From A&S 9.2.28 we require:
0129       //    v < x * eps^1/8
0130       // and from A&S 9.2.29 we require:
0131       //    v^12/10 < 1.5 * x * eps^1/10
0132       // using the former seems to work OK in practice with broadly similar
0133       // error rates either side of the divide for v < 10000.
0134       // At double precision eps^1/8 ~= 0.01.
0135       //
0136       BOOST_MATH_ASSERT(v >= 0);
0137       return (v ? v : 1) < x * 0.004f;
0138 }
0139
0140 template <class T>
0141 inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
0142 {
0143    BOOST_MATH_STD_USING
0144    //
0145    // Determines if x is large enough compared to v to take the asymptotic
0146    // forms above.  From A&S 9.2.28 we require:
0147    //    v < x * eps^1/8
0148    // and from A&S 9.2.29 we require:
0149    //    v^12/10 < 1.5 * x * eps^1/10
0150    // using the former seems to work OK in practice with broadly similar
0151    // error rates either side of the divide for v < 10000.
0152    // At double precision eps^1/8 ~= 0.01.
0153    //
0154    return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
0155 }
0156
0157 template <class T, class Policy>
0158 void temme_asymptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
0159 {
0160    T c = 1;
0161    T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol);
0162    T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
0163    T f = (p - q) / v;
0164    T g_prefix = boost::math::sin_pi(v / 2, pol);
0165    g_prefix *= g_prefix * 2 / v;
0166    T g = f + g_prefix * q;
0167    T h = p;
0168    T c_mult = -x * x / 4;
0169
0170    T y(c * g), y1(c * h);
0171
0172    for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
0173    {
0174       f = (k * f + p + q) / (k*k - v*v);
0175       p /= k - v;
0176       q /= k + v;
0177       c *= c_mult / k;
0178       T c1 = pow(-x * x / 4, T(k)) / factorial<T>(k, pol);
0179       g = f + g_prefix * q;
0180       h = -k * g + p;
0181       y += c * g;
0182       y1 += c * h;
0183       if(c * g / tools::epsilon<T>() < y)
0184          break;
0185    }
0186
0187    *Y = -y;
0188    *Y1 = (-2 / x) * y1;
0189 }
0190
0191 template <class T, class Policy>
0192 T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
0193 {
0194    BOOST_MATH_STD_USING  // ADL of std names
0195    T s = 1;
0196    T mu = 4 * v * v;
0197    T ex = 8 * x;
0198    T num = mu - 1;
0199    T denom = ex;
0200
0201    s -= num / denom;
0202
0203    num *= mu - 9;
0204    denom *= ex * 2;
0205    s += num / denom;
0206
0207    num *= mu - 25;
0208    denom *= ex * 3;
0209    s -= num / denom;
0210
0211    // Try and avoid overflow to the last minute:
0212    T e = exp(x/2);
0213
0214    s = e * (e * s / sqrt(2 * x * constants::pi<T>()));
0215
0216    return (boost::math::isfinite)(s) ?
0217       s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", nullptr, pol);
0218 }
0219
0220 }}} // namespaces
0221
0222 #endif
0223