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
151 if(!need_logs)
152 {
153 gam = boost::math::tgamma(v, pol);
154 p = pow(x / 2, v);
155 if(tools::max_value<T>() * p < gam)
156 {
157 scale /= gam;
158 gam = 1;
159 /*
160 * We can never get here, it would require p < 1/max_value.
161 if(tools::max_value<T>() * p < gam)
162 {
163 return -policies::raise_overflow_error<T>(function, nullptr, pol);
164 }
165 */
166 }
167 prefix = -gam / (constants::pi<T>() * p);
168 }
169 else
170 {
171 gam = boost::math::lgamma(v, pol);
172 p = v * p;
173 prefix = gam - log(constants::pi<T>()) - p;
174 if(tools::log_max_value<T>() < prefix)
175 {
176 prefix -= log(tools::max_value<T>() / 4);
177 scale /= (tools::max_value<T>() / 4);
178 if(tools::log_max_value<T>() < prefix)
179 {
180 return -policies::raise_overflow_error<T>(function, nullptr, pol);
181 }
182 }
183 prefix = -exp(prefix);
184 }
185 bessel_y_small_z_series_term_a<T, Policy> s(v, x);
186 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
187 *pscale = scale;
188
189 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
190
191 policies::check_series_iterations<T>("boost::math::bessel_y_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
192 result *= prefix;
193
194 if(!need_logs)
195 {
196 prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v, pol) * p / constants::pi<T>();
197 }
198 else
199 {
200 int sgn {};
201 prefix = boost::math::lgamma(-v, &sgn, pol) + p;
202 prefix = exp(prefix) * sgn / constants::pi<T>();
203 }
204 bessel_y_small_z_series_term_b<T, Policy> s2(v, x);
205 max_iter = policies::get_max_series_iterations<Policy>();
206
207 T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
208
209 result -= scale * prefix * b;
210 return result;
211}
212
213template <class T, class Policy>
214T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol)
215{
216 //
217 // See http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
218 //
219 // Note that when called we assume that x < epsilon and n is a positive integer.
220 //
221 BOOST_MATH_STD_USING
222 BOOST_MATH_ASSERT(n >= 0);
223 BOOST_MATH_ASSERT((z < policies::get_epsilon<T, Policy>()));
224
225 if(n == 0)
226 {
227 return (2 / constants::pi<T>()) * (log(z / 2) + constants::euler<T>());
228 }
229 else if(n == 1)
230 {
231 return (z / constants::pi<T>()) * log(z / 2)
232 - 2 / (constants::pi<T>() * z)
233 - (z / (2 * constants::pi<T>())) * (1 - 2 * constants::euler<T>());
234 }
235 else if(n == 2)
236 {
237 return (z * z) / (4 * constants::pi<T>()) * log(z / 2)
238 - (4 / (constants::pi<T>() * z * z))
239 - ((z * z) / (8 * constants::pi<T>())) * (T(3)/2 - 2 * constants::euler<T>());
240 }
241 else
242 {
243 #if (defined(__GNUC__) && __GNUC__ == 13)
244 auto p = static_cast<T>(pow(z / 2, T(n)));
245 #else
246 auto p = static_cast<T>(pow(z / 2, n));
247 #endif
248
249 T result = -((boost::math::factorial<T>(n - 1, pol) / constants::pi<T>()));
250 if(p * tools::max_value<T>() < fabs(result))
251 {
252 T div = tools::max_value<T>() / 8;
253 result /= div;
254 *scale /= div;
255 if(p * tools::max_value<T>() < result)
256 {
257 // Impossible to get here??
258 return -policies::raise_overflow_error<T>("bessel_yn_small_z<%1%>(%1%,%1%)", nullptr, pol); // LCOV_EXCL_LINE
259 }
260 }
261 return result / p;
262 }
263}
264
265}}} // namespaces
266
267#endif // BOOST_MATH_BESSEL_JN_SERIES_HPP
268
269

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