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_DERIVATIVES_HPP
7 #define BOOST_MATH_BESSEL_DERIVATIVES_HPP
8 
9 #ifdef _MSC_VER
10 #  pragma once
11 #endif
12 
13 #include <boost/math/special_functions/math_fwd.hpp>
14 #include <boost/math/special_functions/bessel.hpp>
15 #include <boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp>
16 #include <boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp>
17 #include <boost/math/special_functions/detail/bessel_derivatives_linear.hpp>
18 
19 namespace boost{ namespace math{
20 
21 namespace detail{
22 
23 template <class Tag, class T, class Policy>
cyl_bessel_j_prime_imp(T v,T x,const Policy & pol)24 inline T cyl_bessel_j_prime_imp(T v, T x, const Policy& pol)
25 {
26    static const char* const function = "boost::math::cyl_bessel_j_prime<%1%>(%1%,%1%)";
27    BOOST_MATH_STD_USING
28    //
29    // Prevent complex result:
30    //
31    if (x < 0 && floor(v) != v)
32       return boost::math::policies::raise_domain_error<T>(
33          function,
34          "Got x = %1%, but function requires x >= 0", x, pol);
35    //
36    // Special cases for x == 0:
37    //
38    if (x == 0)
39    {
40       if (v == 1)
41          return 0.5;
42       else if (v == -1)
43          return -0.5;
44       else if (floor(v) == v || v > 1)
45          return 0;
46       else return boost::math::policies::raise_domain_error<T>(
47          function,
48          "Got x = %1%, but function is indeterminate for this order", x, pol);
49    }
50    //
51    // Special case for large x: use asymptotic expansion:
52    //
53    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
54       return boost::math::detail::asymptotic_bessel_j_derivative_large_x_2(v, x);
55    //
56    // Special case for small x: use Taylor series:
57    //
58    if ((abs(x) < 5) || (abs(v) > x * x / 4))
59    {
60       bool inversed = false;
61       if (floor(v) == v && v < 0)
62       {
63          v = -v;
64          if (itrunc(v, pol) & 1)
65             inversed = true;
66       }
67       T r = boost::math::detail::bessel_j_derivative_small_z_series(v, x, pol);
68       return inversed ? T(-r) : r;
69    }
70    //
71    // Special case for v == 0:
72    //
73    if (v == 0)
74       return -boost::math::detail::cyl_bessel_j_imp<T>(1, x, Tag(), pol);
75    //
76    // Default case:
77    //
78    return boost::math::detail::bessel_j_derivative_linear(v, x, Tag(), pol);
79 }
80 
81 template <class T, class Policy>
sph_bessel_j_prime_imp(unsigned v,T x,const Policy & pol)82 inline T sph_bessel_j_prime_imp(unsigned v, T x, const Policy& pol)
83 {
84    static const char* const function = "boost::math::sph_bessel_prime<%1%>(%1%,%1%)";
85    //
86    // Prevent complex result:
87    //
88    if (x < 0)
89       return boost::math::policies::raise_domain_error<T>(
90          function,
91          "Got x = %1%, but function requires x >= 0.", x, pol);
92    //
93    // Special case for v == 0:
94    //
95    if (v == 0)
96       return (x == 0) ? boost::math::policies::raise_overflow_error<T>(function, 0, pol)
97          : static_cast<T>(-boost::math::detail::sph_bessel_j_imp<T>(1, x, pol));
98    //
99    // Special case for x == 0 and v > 0:
100    //
101    if (x == 0)
102       return boost::math::policies::raise_domain_error<T>(
103          function,
104          "Got x = %1%, but function is indeterminate for this order", x, pol);
105    //
106    // Default case:
107    //
108    return boost::math::detail::sph_bessel_j_derivative_linear(v, x, pol);
109 }
110 
111 template <class T, class Policy>
cyl_bessel_i_prime_imp(T v,T x,const Policy & pol)112 inline T cyl_bessel_i_prime_imp(T v, T x, const Policy& pol)
113 {
114    static const char* const function = "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)";
115    BOOST_MATH_STD_USING
116    //
117    // Prevent complex result:
118    //
119    if (x < 0 && floor(v) != v)
120       return boost::math::policies::raise_domain_error<T>(
121          function,
122          "Got x = %1%, but function requires x >= 0", x, pol);
123    //
124    // Special cases for x == 0:
125    //
126    if (x == 0)
127    {
128       if (v == 1 || v == -1)
129          return 0.5;
130       else if (floor(v) == v || v > 1)
131          return 0;
132       else return boost::math::policies::raise_domain_error<T>(
133          function,
134          "Got x = %1%, but function is indeterminate for this order", x, pol);
135    }
136    //
137    // Special case for v == 0:
138    //
139    if (v == 0)
140       return boost::math::detail::cyl_bessel_i_imp<T>(1, x, pol);
141    //
142    // Default case:
143    //
144    return boost::math::detail::bessel_i_derivative_linear(v, x, pol);
145 }
146 
147 template <class Tag, class T, class Policy>
cyl_bessel_k_prime_imp(T v,T x,const Policy & pol)148 inline T cyl_bessel_k_prime_imp(T v, T x, const Policy& pol)
149 {
150    //
151    // Prevent complex and indeterminate results:
152    //
153    if (x <= 0)
154       return boost::math::policies::raise_domain_error<T>(
155          "boost::math::cyl_bessel_k_prime<%1%>(%1%,%1%)",
156          "Got x = %1%, but function requires x > 0", x, pol);
157    //
158    // Special case for v == 0:
159    //
160    if (v == 0)
161       return -boost::math::detail::cyl_bessel_k_imp<T>(1, x, Tag(), pol);
162    //
163    // Default case:
164    //
165    return boost::math::detail::bessel_k_derivative_linear(v, x, Tag(), pol);
166 }
167 
168 template <class Tag, class T, class Policy>
cyl_neumann_prime_imp(T v,T x,const Policy & pol)169 inline T cyl_neumann_prime_imp(T v, T x, const Policy& pol)
170 {
171    BOOST_MATH_STD_USING
172    //
173    // Prevent complex and indeterminate results:
174    //
175    if (x <= 0)
176       return boost::math::policies::raise_domain_error<T>(
177          "boost::math::cyl_neumann_prime<%1%>(%1%,%1%)",
178          "Got x = %1%, but function requires x > 0", x, pol);
179    //
180    // Special case for large x: use asymptotic expansion:
181    //
182    if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))
183       return boost::math::detail::asymptotic_bessel_y_derivative_large_x_2(v, x);
184    //
185    // Special case for small x: use Taylor series:
186    //
187    if (v > 0 && floor(v) != v)
188    {
189       const T eps = boost::math::policies::get_epsilon<T, Policy>();
190       if (log(eps / 2) > v * log((x * x) / (v * 4)))
191          return boost::math::detail::bessel_y_derivative_small_z_series(v, x, pol);
192    }
193    //
194    // Special case for v == 0:
195    //
196    if (v == 0)
197       return -boost::math::detail::cyl_neumann_imp<T>(1, x, Tag(), pol);
198    //
199    // Default case:
200    //
201    return boost::math::detail::bessel_y_derivative_linear(v, x, Tag(), pol);
202 }
203 
204 template <class T, class Policy>
sph_neumann_prime_imp(unsigned v,T x,const Policy & pol)205 inline T sph_neumann_prime_imp(unsigned v, T x, const Policy& pol)
206 {
207    //
208    // Prevent complex and indeterminate result:
209    //
210    if (x <= 0)
211       return boost::math::policies::raise_domain_error<T>(
212          "boost::math::sph_neumann_prime<%1%>(%1%,%1%)",
213          "Got x = %1%, but function requires x > 0.", x, pol);
214    //
215    // Special case for v == 0:
216    //
217    if (v == 0)
218       return -boost::math::detail::sph_neumann_imp<T>(1, x, pol);
219    //
220    // Default case:
221    //
222    return boost::math::detail::sph_neumann_derivative_linear(v, x, pol);
223 }
224 
225 } // namespace detail
226 
227 template <class T1, class T2, class Policy>
cyl_bessel_j_prime(T1 v,T2 x,const Policy &)228 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j_prime(T1 v, T2 x, const Policy& /* pol */)
229 {
230    BOOST_FPU_EXCEPTION_GUARD
231    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
232    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
233    typedef typename policies::evaluation<result_type, Policy>::type value_type;
234    typedef typename policies::normalise<
235       Policy,
236       policies::promote_float<false>,
237       policies::promote_double<false>,
238       policies::discrete_quantile<>,
239       policies::assert_undefined<> >::type forwarding_policy;
240    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_j_prime<%1%,%1%>(%1%,%1%)");
241 }
242 
243 template <class T1, class T2>
cyl_bessel_j_prime(T1 v,T2 x)244 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j_prime(T1 v, T2 x)
245 {
246    return cyl_bessel_j_prime(v, x, policies::policy<>());
247 }
248 
249 template <class T, class Policy>
sph_bessel_prime(unsigned v,T x,const Policy &)250 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel_prime(unsigned v, T x, const Policy& /* pol */)
251 {
252    BOOST_FPU_EXCEPTION_GUARD
253    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
254    typedef typename policies::evaluation<result_type, Policy>::type value_type;
255    typedef typename policies::normalise<
256       Policy,
257       policies::promote_float<false>,
258       policies::promote_double<false>,
259       policies::discrete_quantile<>,
260       policies::assert_undefined<> >::type forwarding_policy;
261    return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_bessel_j_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_bessel_j_prime<%1%>(%1%,%1%)");
262 }
263 
264 template <class T>
sph_bessel_prime(unsigned v,T x)265 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel_prime(unsigned v, T x)
266 {
267    return sph_bessel_prime(v, x, policies::policy<>());
268 }
269 
270 template <class T1, class T2, class Policy>
cyl_bessel_i_prime(T1 v,T2 x,const Policy &)271 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i_prime(T1 v, T2 x, const Policy& /* pol */)
272 {
273    BOOST_FPU_EXCEPTION_GUARD
274    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
275    typedef typename policies::evaluation<result_type, Policy>::type value_type;
276    typedef typename policies::normalise<
277       Policy,
278       policies::promote_float<false>,
279       policies::promote_double<false>,
280       policies::discrete_quantile<>,
281       policies::assert_undefined<> >::type forwarding_policy;
282    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_prime_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i_prime<%1%>(%1%,%1%)");
283 }
284 
285 template <class T1, class T2>
cyl_bessel_i_prime(T1 v,T2 x)286 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i_prime(T1 v, T2 x)
287 {
288    return cyl_bessel_i_prime(v, x, policies::policy<>());
289 }
290 
291 template <class T1, class T2, class Policy>
cyl_bessel_k_prime(T1 v,T2 x,const Policy &)292 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k_prime(T1 v, T2 x, const Policy& /* pol */)
293 {
294    BOOST_FPU_EXCEPTION_GUARD
295    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
296    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
297    typedef typename policies::evaluation<result_type, Policy>::type value_type;
298    typedef typename policies::normalise<
299       Policy,
300       policies::promote_float<false>,
301       policies::promote_double<false>,
302       policies::discrete_quantile<>,
303       policies::assert_undefined<> >::type forwarding_policy;
304    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_k_prime<%1%,%1%>(%1%,%1%)");
305 }
306 
307 template <class T1, class T2>
cyl_bessel_k_prime(T1 v,T2 x)308 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k_prime(T1 v, T2 x)
309 {
310    return cyl_bessel_k_prime(v, x, policies::policy<>());
311 }
312 
313 template <class T1, class T2, class Policy>
cyl_neumann_prime(T1 v,T2 x,const Policy &)314 inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann_prime(T1 v, T2 x, const Policy& /* pol */)
315 {
316    BOOST_FPU_EXCEPTION_GUARD
317    typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;
318    typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;
319    typedef typename policies::evaluation<result_type, Policy>::type value_type;
320    typedef typename policies::normalise<
321       Policy,
322       policies::promote_float<false>,
323       policies::promote_double<false>,
324       policies::discrete_quantile<>,
325       policies::assert_undefined<> >::type forwarding_policy;
326    return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_prime_imp<tag_type, value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_neumann_prime<%1%,%1%>(%1%,%1%)");
327 }
328 
329 template <class T1, class T2>
cyl_neumann_prime(T1 v,T2 x)330 inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann_prime(T1 v, T2 x)
331 {
332    return cyl_neumann_prime(v, x, policies::policy<>());
333 }
334 
335 template <class T, class Policy>
sph_neumann_prime(unsigned v,T x,const Policy &)336 inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann_prime(unsigned v, T x, const Policy& /* pol */)
337 {
338    BOOST_FPU_EXCEPTION_GUARD
339    typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;
340    typedef typename policies::evaluation<result_type, Policy>::type value_type;
341    typedef typename policies::normalise<
342       Policy,
343       policies::promote_float<false>,
344       policies::promote_double<false>,
345       policies::discrete_quantile<>,
346       policies::assert_undefined<> >::type forwarding_policy;
347    return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_neumann_prime_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_neumann_prime<%1%>(%1%,%1%)");
348 }
349 
350 template <class T>
sph_neumann_prime(unsigned v,T x)351 inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann_prime(unsigned v, T x)
352 {
353    return sph_neumann_prime(v, x, policies::policy<>());
354 }
355 
356 } // namespace math
357 } // namespace boost
358 
359 #endif // BOOST_MATH_BESSEL_DERIVATIVES_HPP
360