1 // Copyright (c) 2013 Anton Bikineev
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_JY_DERIVATIVES_SERIES_HPP
7 #define BOOST_MATH_BESSEL_JY_DERIVATIVES_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_derivative_small_z_series_term
17 {
18 typedef T result_type;
19
bessel_j_derivative_small_z_series_termboost::math::detail::bessel_j_derivative_small_z_series_term20 bessel_j_derivative_small_z_series_term(T v_, T x)
21 : N(0), v(v_), term(1), mult(x / 2)
22 {
23 mult *= -mult;
24 // iterate if v == 0; otherwise result of
25 // first term is 0 and tools::sum_series stops
26 if (v == 0)
27 iterate();
28 }
operator ()boost::math::detail::bessel_j_derivative_small_z_series_term29 T operator()()
30 {
31 T r = term * (v + 2 * N);
32 iterate();
33 return r;
34 }
35 private:
iterateboost::math::detail::bessel_j_derivative_small_z_series_term36 void iterate()
37 {
38 ++N;
39 term *= mult / (N * (N + v));
40 }
41 unsigned N;
42 T v;
43 T term;
44 T mult;
45 };
46 //
47 // Series evaluation for BesselJ'(v, z) as z -> 0.
48 // It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
49 // Converges rapidly for all z << v.
50 //
51 template <class T, class Policy>
bessel_j_derivative_small_z_series(T v,T x,const Policy & pol)52 inline T bessel_j_derivative_small_z_series(T v, T x, const Policy& pol)
53 {
54 BOOST_MATH_STD_USING
55 T prefix;
56 if (v < boost::math::max_factorial<T>::value)
57 {
58 prefix = pow(x / 2, v - 1) / 2 / boost::math::tgamma(v + 1, pol);
59 }
60 else
61 {
62 prefix = (v - 1) * log(x / 2) - constants::ln_two<T>() - boost::math::lgamma(v + 1, pol);
63 prefix = exp(prefix);
64 }
65 if (0 == prefix)
66 return prefix;
67
68 bessel_j_derivative_small_z_series_term<T, Policy> s(v, x);
69 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
70 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
71 T zero = 0;
72 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
73 #else
74 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
75 #endif
76 boost::math::policies::check_series_iterations<T>("boost::math::bessel_j_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
77 return prefix * result;
78 }
79
80 template <class T, class Policy>
81 struct bessel_y_derivative_small_z_series_term_a
82 {
83 typedef T result_type;
84
bessel_y_derivative_small_z_series_term_aboost::math::detail::bessel_y_derivative_small_z_series_term_a85 bessel_y_derivative_small_z_series_term_a(T v_, T x)
86 : N(0), v(v_)
87 {
88 mult = x / 2;
89 mult *= -mult;
90 term = 1;
91 }
operator ()boost::math::detail::bessel_y_derivative_small_z_series_term_a92 T operator()()
93 {
94 T r = term * (-v + 2 * N);
95 ++N;
96 term *= mult / (N * (N - v));
97 return r;
98 }
99 private:
100 unsigned N;
101 T v;
102 T mult;
103 T term;
104 };
105
106 template <class T, class Policy>
107 struct bessel_y_derivative_small_z_series_term_b
108 {
109 typedef T result_type;
110
bessel_y_derivative_small_z_series_term_bboost::math::detail::bessel_y_derivative_small_z_series_term_b111 bessel_y_derivative_small_z_series_term_b(T v_, T x)
112 : N(0), v(v_)
113 {
114 mult = x / 2;
115 mult *= -mult;
116 term = 1;
117 }
operator ()boost::math::detail::bessel_y_derivative_small_z_series_term_b118 T operator()()
119 {
120 T r = term * (v + 2 * N);
121 ++N;
122 term *= mult / (N * (N + v));
123 return r;
124 }
125 private:
126 unsigned N;
127 T v;
128 T mult;
129 T term;
130 };
131 //
132 // Series form for BesselY' as z -> 0,
133 // It's derivative of 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 //
140 template <class T, class Policy>
bessel_y_derivative_small_z_series(T v,T x,const Policy & pol)141 inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
142 {
143 BOOST_MATH_STD_USING
144 static const char* function = "bessel_y_derivative_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 >= boost::math::max_factorial<T>::value) || (boost::math::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 + 1) * 2;
154 if (boost::math::tools::max_value<T>() * p < gam)
155 {
156 scale /= gam;
157 gam = 1;
158 if (boost::math::tools::max_value<T>() * p < gam)
159 {
160 // This term will overflow to -INF, when combined with the series below it becomes +INF:
161 return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
162 }
163 }
164 prefix = -gam / (boost::math::constants::pi<T>() * p);
165 }
166 else
167 {
168 gam = boost::math::lgamma(v, pol);
169 p = (v + 1) * p + constants::ln_two<T>();
170 prefix = gam - log(boost::math::constants::pi<T>()) - p;
171 if (boost::math::tools::log_max_value<T>() < prefix)
172 {
173 prefix -= log(boost::math::tools::max_value<T>() / 4);
174 scale /= (boost::math::tools::max_value<T>() / 4);
175 if (boost::math::tools::log_max_value<T>() < prefix)
176 {
177 return boost::math::policies::raise_overflow_error<T>(function, 0, pol);
178 }
179 }
180 prefix = -exp(prefix);
181 }
182 bessel_y_derivative_small_z_series_term_a<T, Policy> s(v, x);
183 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
184 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
185 T zero = 0;
186 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
187 #else
188 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
189 #endif
190 boost::math::policies::check_series_iterations<T>("boost::math::bessel_y_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
191 result *= prefix;
192
193 p = pow(x / 2, v - 1) / 2;
194 if (!need_logs)
195 {
196 prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / boost::math::constants::pi<T>();
197 }
198 else
199 {
200 int sgn;
201 prefix = boost::math::lgamma(-v, &sgn, pol) + (v - 1) * log(x / 2) - constants::ln_two<T>();
202 prefix = exp(prefix) * sgn / boost::math::constants::pi<T>();
203 }
204 bessel_y_derivative_small_z_series_term_b<T, Policy> s2(v, x);
205 max_iter = boost::math::policies::get_max_series_iterations<Policy>();
206 #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
207 T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
208 #else
209 T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
210 #endif
211 result += scale * prefix * b;
212 return result;
213 }
214
215 // Calculating of BesselY'(v,x) with small x (x < epsilon) and integer x using derivatives
216 // of formulas in http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
217 // seems to lose precision. Instead using linear combination of regular Bessel is preferred.
218
219 }}} // namespaces
220
221 #endif // BOOST_MATH_BESSEL_JY_DERIVATVIES_SERIES_HPP
222