bessel_prime.hpp 13 KB

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