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) |
5 | |
6 | // |
7 | // This is a partial header, do not include on it's own!!! |
8 | // |
9 | // Contains asymptotic expansions for Bessel J(v,x) and Y(v,x) |
10 | // functions, as x -> INF. |
11 | // |
12 | #ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP |
13 | #define BOOST_MATH_SF_DETAIL_BESSEL_JY_ASYM_HPP |
14 | |
15 | #ifdef _MSC_VER |
16 | #pragma once |
17 | #endif |
18 | |
19 | #include <boost/math/special_functions/factorials.hpp> |
20 | |
21 | namespace boost{ namespace math{ namespace detail{ |
22 | |
23 | template <class T> |
24 | inline T asymptotic_bessel_amplitude(T v, T x) |
25 | { |
26 | // Calculate the amplitude of J(v, x) and Y(v, x) for large |
27 | // x: see A&S 9.2.28. |
28 | BOOST_MATH_STD_USING |
29 | T s = 1; |
30 | T mu = 4 * v * v; |
31 | T txq = 2 * x; |
32 | txq *= txq; |
33 | |
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); |
37 | |
38 | return sqrt(s * 2 / (constants::pi<T>() * x)); |
39 | } |
40 | |
41 | template <class T> |
42 | T asymptotic_bessel_phase_mx(T v, T x) |
43 | { |
44 | // |
45 | // Calculate the phase of J(v, x) and Y(v, x) for large x. |
46 | // See A&S 9.2.29. |
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: |
49 | // |
50 | T mu = 4 * v * v; |
51 | T denom = 4 * x; |
52 | T denom_mult = denom * denom; |
53 | |
54 | T s = 0; |
55 | s += (mu - 1) / (2 * denom); |
56 | denom *= denom_mult; |
57 | s += (mu - 1) * (mu - 25) / (6 * denom); |
58 | denom *= denom_mult; |
59 | s += (mu - 1) * (mu * mu - 114 * mu + 1073) / (5 * denom); |
60 | denom *= denom_mult; |
61 | s += (mu - 1) * (5 * mu * mu * mu - 1535 * mu * mu + 54703 * mu - 375733) / (14 * denom); |
62 | return s; |
63 | } |
64 | |
65 | template <class T> |
66 | inline T asymptotic_bessel_y_large_x_2(T v, T x) |
67 | { |
68 | // See A&S 9.2.19. |
69 | BOOST_MATH_STD_USING |
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); |
75 | // |
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. |
80 | // |
81 | T cx = cos(x); |
82 | T sx = sin(x); |
83 | T ci = cos_pi(v / 2 + 0.25f); |
84 | T si = sin_pi(v / 2 + 0.25f); |
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; |
91 | } |
92 | |
93 | template <class T> |
94 | inline T asymptotic_bessel_j_large_x_2(T v, T x) |
95 | { |
96 | // See A&S 9.2.19. |
97 | BOOST_MATH_STD_USING |
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); |
103 | // |
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. |
108 | // |
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)); |
113 | T cx = cos(x); |
114 | T sx = sin(x); |
115 | T ci = cos_pi(v / 2 + 0.25f); |
116 | T si = sin_pi(v / 2 + 0.25f); |
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; |
120 | } |
121 | |
122 | template <class T> |
123 | inline bool asymptotic_bessel_large_x_limit(int v, const T& x) |
124 | { |
125 | BOOST_MATH_STD_USING |
126 | // |
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: |
129 | // v < x * eps^1/8 |
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. |
135 | // |
136 | BOOST_ASSERT(v >= 0); |
137 | return (v ? v : 1) < x * 0.004f; |
138 | } |
139 | |
140 | template <class T> |
141 | inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x) |
142 | { |
143 | BOOST_MATH_STD_USING |
144 | // |
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: |
147 | // v < x * eps^1/8 |
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. |
153 | // |
154 | return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>()); |
155 | } |
156 | |
157 | template <class T, class Policy> |
158 | void temme_asyptotic_y_small_x(T v, T x, T* Y, T* Y1, const Policy& pol) |
159 | { |
160 | T c = 1; |
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); |
163 | T f = (p - q) / v; |
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; |
167 | T h = p; |
168 | T c_mult = -x * x / 4; |
169 | |
170 | T y(c * g), y1(c * h); |
171 | |
172 | for(int k = 1; k < policies::get_max_series_iterations<Policy>(); ++k) |
173 | { |
174 | f = (k * f + p + q) / (k*k - v*v); |
175 | p /= k - v; |
176 | q /= k + v; |
177 | c *= c_mult / k; |
178 | T c1 = pow(-x * x / 4, k) / factorial<T>(k, pol); |
179 | g = f + g_prefix * q; |
180 | h = -k * g + p; |
181 | y += c * g; |
182 | y1 += c * h; |
183 | if(c * g / tools::epsilon<T>() < y) |
184 | break; |
185 | } |
186 | |
187 | *Y = -y; |
188 | *Y1 = (-2 / x) * y1; |
189 | } |
190 | |
191 | template <class T, class Policy> |
192 | T asymptotic_bessel_i_large_x(T v, T x, const Policy& pol) |
193 | { |
194 | BOOST_MATH_STD_USING // ADL of std names |
195 | T s = 1; |
196 | T mu = 4 * v * v; |
197 | T ex = 8 * x; |
198 | T num = mu - 1; |
199 | T denom = ex; |
200 | |
201 | s -= num / denom; |
202 | |
203 | num *= mu - 9; |
204 | denom *= ex * 2; |
205 | s += num / denom; |
206 | |
207 | num *= mu - 25; |
208 | denom *= ex * 3; |
209 | s -= num / denom; |
210 | |
211 | // Try and avoid overflow to the last minute: |
212 | T e = exp(x/2); |
213 | |
214 | s = e * (e * s / sqrt(2 * x * constants::pi<T>())); |
215 | |
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); |
218 | } |
219 | |
220 | }}} // namespaces |
221 | |
222 | #endif |
223 | |
224 | |