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    if (0 == ampl)
0073       return ampl;
0074    T phase = asymptotic_bessel_phase_mx(v, x);
0075    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
0076    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
0077    //
0078    // Calculate the sine of the phase, using
0079    // sine/cosine addition rules to factor in
0080    // the x - PI(v/2 + 1/4) term not added to the
0081    // phase when we calculated it.
0082    //
0083    T cx = cos(x);
0084    T sx = sin(x);
0085    T ci = boost::math::cos_pi(v / 2 + 0.25f, pol);
0086    T si = boost::math::sin_pi(v / 2 + 0.25f, pol);
0087    T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
0088    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
0089    BOOST_MATH_INSTRUMENT_CODE(cos(x));
0090    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
0091    BOOST_MATH_INSTRUMENT_CODE(sin(x));
0092    return sin_phase * ampl;
0093 }
0094
0095 template <class T, class Policy>
0096 inline T asymptotic_bessel_j_large_x_2(T v, T x, const Policy& pol)
0097 {
0098    // See A&S 9.2.19.
0099    BOOST_MATH_STD_USING
0100    // Get the phase and amplitude:
0101    T ampl = asymptotic_bessel_amplitude(v, x);
0102    if (0 == ampl)
0103       return ampl;  // shortcut.
0104    T phase = asymptotic_bessel_phase_mx(v, x);
0105    BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
0106    BOOST_MATH_INSTRUMENT_VARIABLE(phase);
0107    //
0108    // Calculate the sine of the phase, using
0109    // sine/cosine addition rules to factor in
0110    // the x - PI(v/2 + 1/4) term not added to the
0111    // phase when we calculated it.
0112    //
0113    BOOST_MATH_INSTRUMENT_CODE(cos(phase));
0114    BOOST_MATH_INSTRUMENT_CODE(cos(x));
0115    BOOST_MATH_INSTRUMENT_CODE(sin(phase));
0116    BOOST_MATH_INSTRUMENT_CODE(sin(x));
0117    T cx = cos(x);
0118    T sx = sin(x);
0119    T ci = boost::math::cos_pi(v / 2 + 0.25f, pol);
0120    T si = boost::math::sin_pi(v / 2 + 0.25f, pol);
0121    T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
0122    BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
0123    return sin_phase * ampl;
0124 }
0125
0126 template <class T>
0127 inline bool asymptotic_bessel_large_x_limit(int v, const T& x)
0128 {
0129    BOOST_MATH_STD_USING
0130       //
0131       // Determines if x is large enough compared to v to take the asymptotic
0132       // forms above.  From A&S 9.2.28 we require:
0133       //    v < x * eps^1/8
0134       // and from A&S 9.2.29 we require:
0135       //    v^12/10 < 1.5 * x * eps^1/10
0136       // using the former seems to work OK in practice with broadly similar
0137       // error rates either side of the divide for v < 10000.
0138       // At double precision eps^1/8 ~= 0.01.
0139       //
0140       BOOST_MATH_ASSERT(v >= 0);
0141       return (v ? v : 1) < x * 0.004f;
0142 }
0143
0144 template <class T>
0145 inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
0146 {
0147    BOOST_MATH_STD_USING
0148    //
0149    // Determines if x is large enough compared to v to take the asymptotic
0150    // forms above.  From A&S 9.2.28 we require:
0151    //    v < x * eps^1/8
0152    // and from A&S 9.2.29 we require:
0153    //    v^12/10 < 1.5 * x * eps^1/10
0154    // using the former seems to work OK in practice with broadly similar
0155    // error rates either side of the divide for v < 10000.
0156    // At double precision eps^1/8 ~= 0.01.
0157    //
0158    return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
0159 }
0160
0161 template <class T, class Policy>
0162 void temme_asymptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
0163 {
0164    T c = 1;
0165    T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol);
0166    T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
0167    T f = (p - q) / v;
0168    T g_prefix = boost::math::sin_pi(v / 2, pol);
0169    g_prefix *= g_prefix * 2 / v;
0170    T g = f + g_prefix * q;
0171    T h = p;
0172    T c_mult = -x * x / 4;
0173
0174    T y(c * g), y1(c * h);
0175
0176    for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
0177    {
0178       f = (k * f + p + q) / (k*k - v*v);
0179       p /= k - v;
0180       q /= k + v;
0181       c *= c_mult / k;
0182       T c1 = pow(-x * x / 4, T(k)) / factorial<T>(k, pol);
0183       g = f + g_prefix * q;
0184       h = -k * g + p;
0185       y += c * g;
0186       y1 += c * h;
0187       if(c * g / tools::epsilon<T>() < y)
0188          break;
0189    }
0190
0191    *Y = -y;
0192    *Y1 = (-2 / x) * y1;
0193 }
0194
0195 template <class T, class Policy>
0196 T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
0197 {
0198    BOOST_MATH_STD_USING  // ADL of std names
0199    T s = 1;
0200    T mu = 4 * v * v;
0201    T ex = 8 * x;
0202    T num = mu - 1;
0203    T denom = ex;
0204
0205    s -= num / denom;
0206
0207    num *= mu - 9;
0208    denom *= ex * 2;
0209    s += num / denom;
0210
0211    num *= mu - 25;
0212    denom *= ex * 3;
0213    s -= num / denom;
0214
0215    // Try and avoid overflow to the last minute:
0216    T e = exp(x/2);
0217
0218    s = e * (e * s / sqrt(2 * x * constants::pi<T>()));
0219
0220    return (boost::math::isfinite)(s) ?
0221       s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", nullptr, pol);
0222 }
0223
0224 }}} // namespaces
0225
0226 #endif
0227