]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright (c) 2011 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 | #ifndef BOOST_MATH_BESSEL_JN_SERIES_HPP | |
7 | #define BOOST_MATH_BESSEL_JN_SERIES_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | namespace boost { namespace math { namespace detail{ | |
14 | ||
15 | template <class T, class Policy> | |
16 | struct bessel_j_small_z_series_term | |
17 | { | |
18 | typedef T result_type; | |
19 | ||
20 | bessel_j_small_z_series_term(T v_, T x) | |
21 | : N(0), v(v_) | |
22 | { | |
23 | BOOST_MATH_STD_USING | |
24 | mult = x / 2; | |
25 | mult *= -mult; | |
26 | term = 1; | |
27 | } | |
28 | T operator()() | |
29 | { | |
30 | T r = term; | |
31 | ++N; | |
32 | term *= mult / (N * (N + v)); | |
33 | return r; | |
34 | } | |
35 | private: | |
36 | unsigned N; | |
37 | T v; | |
38 | T mult; | |
39 | T term; | |
40 | }; | |
41 | // | |
42 | // Series evaluation for BesselJ(v, z) as z -> 0. | |
43 | // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/ | |
44 | // Converges rapidly for all z << v. | |
45 | // | |
46 | template <class T, class Policy> | |
47 | inline T bessel_j_small_z_series(T v, T x, const Policy& pol) | |
48 | { | |
49 | BOOST_MATH_STD_USING | |
50 | T prefix; | |
51 | if(v < max_factorial<T>::value) | |
52 | { | |
53 | prefix = pow(x / 2, v) / boost::math::tgamma(v+1, pol); | |
54 | } | |
55 | else | |
56 | { | |
57 | prefix = v * log(x / 2) - boost::math::lgamma(v+1, pol); | |
58 | prefix = exp(prefix); | |
59 | } | |
60 | if(0 == prefix) | |
61 | return prefix; | |
62 | ||
63 | bessel_j_small_z_series_term<T, Policy> s(v, x); | |
64 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
65 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
66 | T zero = 0; | |
67 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
68 | #else | |
69 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
70 | #endif | |
71 | policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol); | |
72 | return prefix * result; | |
73 | } | |
74 | ||
75 | template <class T, class Policy> | |
76 | struct bessel_y_small_z_series_term_a | |
77 | { | |
78 | typedef T result_type; | |
79 | ||
80 | bessel_y_small_z_series_term_a(T v_, T x) | |
81 | : N(0), v(v_) | |
82 | { | |
83 | BOOST_MATH_STD_USING | |
84 | mult = x / 2; | |
85 | mult *= -mult; | |
86 | term = 1; | |
87 | } | |
88 | T operator()() | |
89 | { | |
90 | BOOST_MATH_STD_USING | |
91 | T r = term; | |
92 | ++N; | |
93 | term *= mult / (N * (N - v)); | |
94 | return r; | |
95 | } | |
96 | private: | |
97 | unsigned N; | |
98 | T v; | |
99 | T mult; | |
100 | T term; | |
101 | }; | |
102 | ||
103 | template <class T, class Policy> | |
104 | struct bessel_y_small_z_series_term_b | |
105 | { | |
106 | typedef T result_type; | |
107 | ||
108 | bessel_y_small_z_series_term_b(T v_, T x) | |
109 | : N(0), v(v_) | |
110 | { | |
111 | BOOST_MATH_STD_USING | |
112 | mult = x / 2; | |
113 | mult *= -mult; | |
114 | term = 1; | |
115 | } | |
116 | T operator()() | |
117 | { | |
118 | T r = term; | |
119 | ++N; | |
120 | term *= mult / (N * (N + v)); | |
121 | return r; | |
122 | } | |
123 | private: | |
124 | unsigned N; | |
125 | T v; | |
126 | T mult; | |
127 | T term; | |
128 | }; | |
129 | // | |
130 | // Series form for BesselY as z -> 0, | |
131 | // see: http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/01/0003/ | |
132 | // This series is only useful when the second term is small compared to the first | |
133 | // otherwise we get catestrophic cancellation errors. | |
134 | // | |
135 | // Approximating tgamma(v) by v^v, and assuming |tgamma(-z)| < eps we end up requiring: | |
136 | // eps/2 * v^v(x/2)^-v > (x/2)^v or log(eps/2) > v log((x/2)^2/v) | |
137 | // | |
138 | template <class T, class Policy> | |
139 | inline T bessel_y_small_z_series(T v, T x, T* pscale, const Policy& pol) | |
140 | { | |
141 | BOOST_MATH_STD_USING | |
142 | static const char* function = "bessel_y_small_z_series<%1%>(%1%,%1%)"; | |
143 | T prefix; | |
144 | T gam; | |
145 | T p = log(x / 2); | |
146 | T scale = 1; | |
147 | bool need_logs = (v >= max_factorial<T>::value) || (tools::log_max_value<T>() / v < fabs(p)); | |
148 | if(!need_logs) | |
149 | { | |
150 | gam = boost::math::tgamma(v, pol); | |
151 | p = pow(x / 2, v); | |
152 | if(tools::max_value<T>() * p < gam) | |
153 | { | |
154 | scale /= gam; | |
155 | gam = 1; | |
156 | if(tools::max_value<T>() * p < gam) | |
157 | { | |
158 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
159 | } | |
160 | } | |
161 | prefix = -gam / (constants::pi<T>() * p); | |
162 | } | |
163 | else | |
164 | { | |
165 | gam = boost::math::lgamma(v, pol); | |
166 | p = v * p; | |
167 | prefix = gam - log(constants::pi<T>()) - p; | |
168 | if(tools::log_max_value<T>() < prefix) | |
169 | { | |
170 | prefix -= log(tools::max_value<T>() / 4); | |
171 | scale /= (tools::max_value<T>() / 4); | |
172 | if(tools::log_max_value<T>() < prefix) | |
173 | { | |
174 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
175 | } | |
176 | } | |
177 | prefix = -exp(prefix); | |
178 | } | |
179 | bessel_y_small_z_series_term_a<T, Policy> s(v, x); | |
180 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
181 | *pscale = scale; | |
182 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
183 | T zero = 0; | |
184 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
185 | #else | |
186 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
187 | #endif | |
188 | policies::check_series_iterations<T>("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol); | |
189 | result *= prefix; | |
190 | ||
191 | if(!need_logs) | |
192 | { | |
193 | prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / constants::pi<T>(); | |
194 | } | |
195 | else | |
196 | { | |
197 | int sgn; | |
198 | prefix = boost::math::lgamma(-v, &sgn, pol) + p; | |
199 | prefix = exp(prefix) * sgn / constants::pi<T>(); | |
200 | } | |
201 | bessel_y_small_z_series_term_b<T, Policy> s2(v, x); | |
202 | max_iter = policies::get_max_series_iterations<Policy>(); | |
203 | #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) | |
204 | T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero); | |
205 | #else | |
206 | T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter); | |
207 | #endif | |
208 | result -= scale * prefix * b; | |
209 | return result; | |
210 | } | |
211 | ||
212 | template <class T, class Policy> | |
213 | T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol) | |
214 | { | |
215 | // | |
216 | // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/ | |
217 | // | |
218 | // Note that when called we assume that x < epsilon and n is a positive integer. | |
219 | // | |
220 | BOOST_MATH_STD_USING | |
221 | BOOST_ASSERT(n >= 0); | |
222 | BOOST_ASSERT((z < policies::get_epsilon<T, Policy>())); | |
223 | ||
224 | if(n == 0) | |
225 | { | |
226 | return (2 / constants::pi<T>()) * (log(z / 2) + constants::euler<T>()); | |
227 | } | |
228 | else if(n == 1) | |
229 | { | |
230 | return (z / constants::pi<T>()) * log(z / 2) | |
231 | - 2 / (constants::pi<T>() * z) | |
232 | - (z / (2 * constants::pi<T>())) * (1 - 2 * constants::euler<T>()); | |
233 | } | |
234 | else if(n == 2) | |
235 | { | |
236 | return (z * z) / (4 * constants::pi<T>()) * log(z / 2) | |
237 | - (4 / (constants::pi<T>() * z * z)) | |
238 | - ((z * z) / (8 * constants::pi<T>())) * (T(3)/2 - 2 * constants::euler<T>()); | |
239 | } | |
240 | else | |
241 | { | |
242 | T p = pow(z / 2, n); | |
243 | T result = -((boost::math::factorial<T>(n - 1) / constants::pi<T>())); | |
244 | if(p * tools::max_value<T>() < result) | |
245 | { | |
246 | T div = tools::max_value<T>() / 8; | |
247 | result /= div; | |
248 | *scale /= div; | |
249 | if(p * tools::max_value<T>() < result) | |
250 | { | |
251 | return -policies::raise_overflow_error<T>("bessel_yn_small_z<%1%>(%1%,%1%)", 0, pol); | |
252 | } | |
253 | } | |
254 | return result / p; | |
255 | } | |
256 | } | |
257 | ||
258 | }}} // namespaces | |
259 | ||
260 | #endif // BOOST_MATH_BESSEL_JN_SERIES_HPP | |
261 |