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)
6 #ifndef BOOST_MATH_BESSEL_JN_SERIES_HPP
7 #define BOOST_MATH_BESSEL_JN_SERIES_HPP
13 namespace boost { namespace math { namespace detail{
15 template <class T, class Policy>
16 struct bessel_j_small_z_series_term
18 typedef T result_type;
20 bessel_j_small_z_series_term(T v_, T x)
32 term *= mult / (N * (N + v));
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.
46 template <class T, class Policy>
47 inline T bessel_j_small_z_series(T v, T x, const Policy& pol)
51 if(v < max_factorial<T>::value)
53 prefix = pow(x / 2, v) / boost::math::tgamma(v+1, pol);
57 prefix = v * log(x / 2) - boost::math::lgamma(v+1, pol);
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))
67 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
69 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
71 policies::check_series_iterations<T>("boost::math::bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
72 return prefix * result;
75 template <class T, class Policy>
76 struct bessel_y_small_z_series_term_a
78 typedef T result_type;
80 bessel_y_small_z_series_term_a(T v_, T x)
93 term *= mult / (N * (N - v));
103 template <class T, class Policy>
104 struct bessel_y_small_z_series_term_b
106 typedef T result_type;
108 bessel_y_small_z_series_term_b(T v_, T x)
120 term *= mult / (N * (N + v));
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.
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)
138 template <class T, class Policy>
139 inline T bessel_y_small_z_series(T v, T x, T* pscale, const Policy& pol)
142 static const char* function = "bessel_y_small_z_series<%1%>(%1%,%1%)";
147 bool need_logs = (v >= max_factorial<T>::value) || (tools::log_max_value<T>() / v < fabs(p));
150 gam = boost::math::tgamma(v, pol);
152 if(tools::max_value<T>() * p < gam)
156 if(tools::max_value<T>() * p < gam)
158 return -policies::raise_overflow_error<T>(function, 0, pol);
161 prefix = -gam / (constants::pi<T>() * p);
165 gam = boost::math::lgamma(v, pol);
167 prefix = gam - log(constants::pi<T>()) - p;
168 if(tools::log_max_value<T>() < prefix)
170 prefix -= log(tools::max_value<T>() / 4);
171 scale /= (tools::max_value<T>() / 4);
172 if(tools::log_max_value<T>() < prefix)
174 return -policies::raise_overflow_error<T>(function, 0, pol);
177 prefix = -exp(prefix);
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>();
182 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
184 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
186 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
188 policies::check_series_iterations<T>("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
193 prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / constants::pi<T>();
198 prefix = boost::math::lgamma(-v, &sgn, pol) + p;
199 prefix = exp(prefix) * sgn / constants::pi<T>();
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);
206 T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
208 result -= scale * prefix * b;
212 template <class T, class Policy>
213 T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol)
216 // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
218 // Note that when called we assume that x < epsilon and n is a positive integer.
221 BOOST_ASSERT(n >= 0);
222 BOOST_ASSERT((z < policies::get_epsilon<T, Policy>()));
226 return (2 / constants::pi<T>()) * (log(z / 2) + constants::euler<T>());
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>());
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>());
243 T result = -((boost::math::factorial<T>(n - 1) / constants::pi<T>()));
244 if(p * tools::max_value<T>() < result)
246 T div = tools::max_value<T>() / 8;
249 if(p * tools::max_value<T>() < result)
251 return -policies::raise_overflow_error<T>("bessel_yn_small_z<%1%>(%1%,%1%)", 0, pol);
260 #endif // BOOST_MATH_BESSEL_JN_SERIES_HPP