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
13namespace boost { namespace math { namespace detail{
14
15template <class T, class Policy>
16struct 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 }
35private:
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//
46template <class T, class Policy>
47inline 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
75template <class T, class Policy>
76struct 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 }
96private:
97 unsigned N;
98 T v;
99 T mult;
100 T term;
101};
102
103template <class T, class Policy>
104struct 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 }
123private:
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 catastrophic 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//
138template <class T, class Policy>
139inline 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
212template <class T, class Policy>
213T 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
262

source code of include/boost/math/special_functions/detail/bessel_jy_series.hpp