1 // Copyright (c) 2007 John Maddock
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
7 // This is a partial header, do not include on it's own!!!
9 // Contains asymptotic expansions for Bessel J(v,x) and Y(v,x)
10 // functions, as x -> INF.
12 #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
13 #define BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP
19 #include <boost/math/special_functions/factorials.hpp>
21 namespace boost{ namespace math{ namespace detail{
24 inline T asymptotic_bessel_amplitude(T v, T x)
26 // Calculate the amplitude of J(v, x) and Y(v, x) for large
34 s += (mu - 1) / (2 * txq);
35 s += 3 * (mu - 1) * (mu - 9) / (txq * txq * 8);
36 s += 15 * (mu - 1) * (mu - 9) * (mu - 25) / (txq * txq * txq * 8 * 6);
38 return sqrt(s * 2 / (constants::pi<T>() * x));
42 T asymptotic_bessel_phase_mx(T v, T x)
45 // Calculate the phase of J(v, x) and Y(v, x) for large x.
47 // Note that the result returned is the phase less (x - PI(v/2 + 1/4))
48 // which we'll factor in later when we calculate the sines/cosines of the result:
52 T denom_mult = denom * denom;
55 s += (mu - 1) / (2 * denom);
57 s += (mu - 1) * (mu - 25) / (6 * denom);
59 s += (mu - 1) * (mu * mu - 114 * mu + 1073) / (5 * denom);
61 s += (mu - 1) * (5 * mu * mu * mu - 1535 * mu * mu + 54703 * mu - 375733) / (14 * denom);
65 template <class T, class Policy>
66 inline T asymptotic_bessel_y_large_x_2(T v, T x, const Policy& pol)
70 // Get the phase and amplitude:
71 T ampl = asymptotic_bessel_amplitude(v, x);
72 T phase = asymptotic_bessel_phase_mx(v, x);
73 BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
74 BOOST_MATH_INSTRUMENT_VARIABLE(phase);
76 // Calculate the sine of the phase, using
77 // sine/cosine addition rules to factor in
78 // the x - PI(v/2 + 1/4) term not added to the
79 // phase when we calculated it.
83 T ci = boost::math::cos_pi(v / 2 + 0.25f, pol);
84 T si = boost::math::sin_pi(v / 2 + 0.25f, pol);
85 T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
86 BOOST_MATH_INSTRUMENT_CODE(sin(phase));
87 BOOST_MATH_INSTRUMENT_CODE(cos(x));
88 BOOST_MATH_INSTRUMENT_CODE(cos(phase));
89 BOOST_MATH_INSTRUMENT_CODE(sin(x));
90 return sin_phase * ampl;
93 template <class T, class Policy>
94 inline T asymptotic_bessel_j_large_x_2(T v, T x, const Policy& pol)
98 // Get the phase and amplitude:
99 T ampl = asymptotic_bessel_amplitude(v, x);
100 T phase = asymptotic_bessel_phase_mx(v, x);
101 BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
102 BOOST_MATH_INSTRUMENT_VARIABLE(phase);
104 // Calculate the sine of the phase, using
105 // sine/cosine addition rules to factor in
106 // the x - PI(v/2 + 1/4) term not added to the
107 // phase when we calculated it.
109 BOOST_MATH_INSTRUMENT_CODE(cos(phase));
110 BOOST_MATH_INSTRUMENT_CODE(cos(x));
111 BOOST_MATH_INSTRUMENT_CODE(sin(phase));
112 BOOST_MATH_INSTRUMENT_CODE(sin(x));
115 T ci = boost::math::cos_pi(v / 2 + 0.25f, pol);
116 T si = boost::math::sin_pi(v / 2 + 0.25f, pol);
117 T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
118 BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
119 return sin_phase * ampl;
123 inline bool asymptotic_bessel_large_x_limit(int v, const T& x)
127 // Determines if x is large enough compared to v to take the asymptotic
128 // forms above. From A&S 9.2.28 we require:
130 // and from A&S 9.2.29 we require:
131 // v^12/10 < 1.5 * x * eps^1/10
132 // using the former seems to work OK in practice with broadly similar
133 // error rates either side of the divide for v < 10000.
134 // At double precision eps^1/8 ~= 0.01.
136 BOOST_MATH_ASSERT(v >= 0);
137 return (v ? v : 1) < x * 0.004f;
141 inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
145 // Determines if x is large enough compared to v to take the asymptotic
146 // forms above. From A&S 9.2.28 we require:
148 // and from A&S 9.2.29 we require:
149 // v^12/10 < 1.5 * x * eps^1/10
150 // using the former seems to work OK in practice with broadly similar
151 // error rates either side of the divide for v < 10000.
152 // At double precision eps^1/8 ~= 0.01.
154 return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
157 template <class T, class Policy>
158 void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol)
161 T p = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, -v) / boost::math::tgamma(1 - v, pol);
162 T q = (v / boost::math::sin_pi(v, pol)) * pow(x / 2, v) / boost::math::tgamma(1 + v, pol);
164 T g_prefix = boost::math::sin_pi(v / 2, pol);
165 g_prefix *= g_prefix * 2 / v;
166 T g = f + g_prefix * q;
168 T c_mult = -x * x / 4;
170 T y(c * g), y1(c * h);
172 for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k)
174 f = (k * f + p + q) / (k*k - v*v);
178 T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol);
179 g = f + g_prefix * q;
183 if(c * g / tools::epsilon<T>() < y)
191 template <class T, class Policy>
192 T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol)
194 BOOST_MATH_STD_USING // ADL of std names
211 // Try and avoid overflow to the last minute:
214 s = e * (e * s / sqrt(2 * x * constants::pi<T>()));
216 return (boost::math::isfinite)(s) ?
217 s : policies::raise_overflow_error<T>("boost::math::asymptotic_bessel_i_large_x<%1%>(%1%,%1%)", 0, pol);