| // Copyright (c) 2007 John Maddock |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. (See accompanying file |
| // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| // |
| // This is a partial header, do not include on it's own!!! |
| // |
| // Contains asymptotic expansions for Bessel J(v,x) and Y(v,x) |
| // functions, as x -> INF. |
| // |
| #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP |
| #define BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #include <boost/math/special_functions/factorials.hpp> |
| |
| namespace boost{ namespace math{ namespace detail{ |
| |
| template <class T> |
| inline T asymptotic_bessel_j_large_x_P(T v, T x) |
| { |
| // A&S 9.2.9 |
| T s = 1; |
| T mu = 4 * v * v; |
| T ez2 = 8 * x; |
| ez2 *= ez2; |
| s -= (mu-1) * (mu-9) / (2 * ez2); |
| s += (mu-1) * (mu-9) * (mu-25) * (mu - 49) / (24 * ez2 * ez2); |
| return s; |
| } |
| |
| template <class T> |
| inline T asymptotic_bessel_j_large_x_Q(T v, T x) |
| { |
| // A&S 9.2.10 |
| T s = 0; |
| T mu = 4 * v * v; |
| T ez = 8*x; |
| s += (mu-1) / ez; |
| s -= (mu-1) * (mu-9) * (mu-25) / (6 * ez*ez*ez); |
| return s; |
| } |
| |
| template <class T> |
| inline T asymptotic_bessel_j_large_x(T v, T x) |
| { |
| // |
| // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/ |
| // |
| // Also A&S 9.2.5 |
| // |
| BOOST_MATH_STD_USING // ADL of std names |
| T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4; |
| return sqrt(2 / (constants::pi<T>() * x)) |
| * (asymptotic_bessel_j_large_x_P(v, x) * cos(chi) |
| - asymptotic_bessel_j_large_x_Q(v, x) * sin(chi)); |
| } |
| |
| template <class T> |
| inline T asymptotic_bessel_y_large_x(T v, T x) |
| { |
| // |
| // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/ |
| // |
| // Also A&S 9.2.5 |
| // |
| BOOST_MATH_STD_USING // ADL of std names |
| T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4; |
| return sqrt(2 / (constants::pi<T>() * x)) |
| * (asymptotic_bessel_j_large_x_P(v, x) * sin(chi) |
| - asymptotic_bessel_j_large_x_Q(v, x) * cos(chi)); |
| } |
| |
| template <class T> |
| inline T asymptotic_bessel_amplitude(T v, T x) |
| { |
| // Calculate the amplitude of J(v, x) and Y(v, x) for large |
| // x: see A&S 9.2.28. |
| BOOST_MATH_STD_USING |
| T s = 1; |
| T mu = 4 * v * v; |
| T txq = 2 * x; |
| txq *= txq; |
| |
| s += (mu - 1) / (2 * txq); |
| s += 3 * (mu - 1) * (mu - 9) / (txq * txq * 8); |
| s += 15 * (mu - 1) * (mu - 9) * (mu - 25) / (txq * txq * txq * 8 * 6); |
| |
| return sqrt(s * 2 / (constants::pi<T>() * x)); |
| } |
| |
| template <class T> |
| T asymptotic_bessel_phase_mx(T v, T x) |
| { |
| // |
| // Calculate the phase of J(v, x) and Y(v, x) for large x. |
| // See A&S 9.2.29. |
| // Note that the result returned is the phase less x. |
| // |
| T mu = 4 * v * v; |
| T denom = 4 * x; |
| T denom_mult = denom * denom; |
| |
| T s = -constants::pi<T>() * (v / 2 + 0.25f); |
| s += (mu - 1) / (2 * denom); |
| denom *= denom_mult; |
| s += (mu - 1) * (mu - 25) / (6 * denom); |
| denom *= denom_mult; |
| s += (mu - 1) * (mu * mu - 114 * mu + 1073) / (5 * denom); |
| denom *= denom_mult; |
| s += (mu - 1) * (5 * mu * mu * mu - 1535 * mu * mu + 54703 * mu - 375733) / (14 * denom); |
| return s; |
| } |
| |
| template <class T> |
| inline T asymptotic_bessel_y_large_x_2(T v, T x) |
| { |
| // See A&S 9.2.19. |
| BOOST_MATH_STD_USING |
| // Get the phase and amplitude: |
| T ampl = asymptotic_bessel_amplitude(v, x); |
| T phase = asymptotic_bessel_phase_mx(v, x); |
| // |
| // Calculate the sine of the phase, using: |
| // sin(x+p) = sin(x)cos(p) + cos(x)sin(p) |
| // |
| T sin_phase = sin(phase) * cos(x) + cos(phase) * sin(x); |
| return sin_phase * ampl; |
| } |
| |
| template <class T> |
| inline T asymptotic_bessel_j_large_x_2(T v, T x) |
| { |
| // See A&S 9.2.19. |
| BOOST_MATH_STD_USING |
| // Get the phase and amplitude: |
| T ampl = asymptotic_bessel_amplitude(v, x); |
| T phase = asymptotic_bessel_phase_mx(v, x); |
| // |
| // Calculate the sine of the phase, using: |
| // cos(x+p) = cos(x)cos(p) - sin(x)sin(p) |
| // |
| T sin_phase = cos(phase) * cos(x) - sin(phase) * sin(x); |
| return sin_phase * ampl; |
| } |
| |
| // |
| // Various limits for the J and Y asymptotics |
| // (the asympotic expansions are safe to use if |
| // x is less than the limit given). |
| // We assume that if we don't use these expansions then the |
| // error will likely be >100eps, so the limits given are chosen |
| // to lead to < 100eps truncation error. |
| // |
| template <class T> |
| inline T asymptotic_bessel_y_limit(const mpl::int_<0>&) |
| { |
| // default case: |
| BOOST_MATH_STD_USING |
| return 2.25 / pow(100 * tools::epsilon<T>() / T(0.001f), T(0.2f)); |
| } |
| template <class T> |
| inline T asymptotic_bessel_y_limit(const mpl::int_<53>&) |
| { |
| // double case: |
| return 304 /*780*/; |
| } |
| template <class T> |
| inline T asymptotic_bessel_y_limit(const mpl::int_<64>&) |
| { |
| // 80-bit extended-double case: |
| return 1552 /*3500*/; |
| } |
| template <class T> |
| inline T asymptotic_bessel_y_limit(const mpl::int_<113>&) |
| { |
| // 128-bit long double case: |
| return 1245243 /*3128000*/; |
| } |
| |
| template <class T, class Policy> |
| struct bessel_asymptotic_tag |
| { |
| typedef typename policies::precision<T, Policy>::type precision_type; |
| typedef typename mpl::if_< |
| mpl::or_< |
| mpl::equal_to<precision_type, mpl::int_<0> >, |
| mpl::greater<precision_type, mpl::int_<113> > >, |
| mpl::int_<0>, |
| typename mpl::if_< |
| mpl::greater<precision_type, mpl::int_<64> >, |
| mpl::int_<113>, |
| typename mpl::if_< |
| mpl::greater<precision_type, mpl::int_<53> >, |
| mpl::int_<64>, |
| mpl::int_<53> |
| >::type |
| >::type |
| >::type type; |
| }; |
| |
| template <class T> |
| inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<0>&) |
| { |
| // default case: |
| BOOST_MATH_STD_USING |
| T v2 = (std::max)(T(3), T(v * v)); |
| return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f)); |
| } |
| template <class T> |
| inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<53>&) |
| { |
| // double case: |
| T v2 = (std::max)(T(3), v * v); |
| return v2 * 33 /*73*/; |
| } |
| template <class T> |
| inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<64>&) |
| { |
| // 80-bit extended-double case: |
| T v2 = (std::max)(T(3), v * v); |
| return v2 * 121 /*266*/; |
| } |
| template <class T> |
| inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<113>&) |
| { |
| // 128-bit long double case: |
| T v2 = (std::max)(T(3), v * v); |
| return v2 * 39154 /*85700*/; |
| } |
| |
| template <class T, class Policy> |
| void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol) |
| { |
| T c = 1; |
| T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol); |
| T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol); |
| T f = (p - q) / v; |
| T g_prefix = boost::math::sin_pi(v / 2, pol); |
| g_prefix *= g_prefix * 2 / v; |
| T g = f + g_prefix * q; |
| T h = p; |
| T c_mult = -x * x / 4; |
| |
| T y(c * g), y1(c * h); |
| |
| for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k) |
| { |
| f = (k * f + p + q) / (k*k - v*v); |
| p /= k - v; |
| q /= k + v; |
| c *= c_mult / k; |
| T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol); |
| g = f + g_prefix * q; |
| h = -k * g + p; |
| y += c * g; |
| y1 += c * h; |
| if(c * g / tools::epsilon<T>() < y) |
| break; |
| } |
| |
| *Y = -y; |
| *Y1 = (-2 / x) * y1; |
| } |
| |
| template <class T, class Policy> |
| T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol) |
| { |
| BOOST_MATH_STD_USING // ADL of std names |
| T s = 1; |
| T mu = 4 * v * v; |
| T ex = 8 * x; |
| T num = mu - 1; |
| T denom = ex; |
| |
| s -= num / denom; |
| |
| num *= mu - 9; |
| denom *= ex * 2; |
| s += num / denom; |
| |
| num *= mu - 25; |
| denom *= ex * 3; |
| s -= num / denom; |
| |
| // Try and avoid overflow to the last minute: |
| T e = exp(x/2); |
| |
| s = e * (e * s / sqrt(2 * x * constants::pi<T>())); |
| |
| return (boost::math::isfinite)(s) ? |
| s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol); |
| } |
| |
| }}} // namespaces |
| |
| #endif |
| |