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

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