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