344 lines · plain
1// Copyright (c) 2013 Anton Bikineev2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5//6#ifndef BOOST_MATH_BESSEL_DERIVATIVES_HPP7#define BOOST_MATH_BESSEL_DERIVATIVES_HPP8 9#ifdef _MSC_VER10# pragma once11#endif12 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 19namespace boost{ namespace math{20 21namespace detail{22 23template <class Tag, class T, class Policy>24inline 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_USING28 //29 // Prevent complex result:30 //31 if ((x < 0) && (floor(v) != v))32 return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x >= 0", x, pol);33 //34 // Special cases for x == 0:35 //36 if (x == 0)37 {38 if (v == 1)39 return static_cast<T>(0.5);40 else if (v == -1)41 return static_cast<T>(-0.5);42 else if (floor(v) == v || v > 1)43 return 0;44 else45 return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function is indeterminate for this order", x, pol);46 }47 //48 // Special case for large x: use asymptotic expansion:49 //50 if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))51 return boost::math::detail::asymptotic_bessel_j_derivative_large_x_2(v, x, pol);52 //53 // Special case for small x: use Taylor series:54 //55 if ((abs(x) < 5) || (abs(v) > x * x / 4))56 {57 bool inversed = false;58 if (floor(v) == v && v < 0)59 {60 v = -v;61 if (itrunc(v, pol) & 1)62 inversed = true;63 }64 T r = boost::math::detail::bessel_j_derivative_small_z_series(v, x, pol);65 return inversed ? T(-r) : r;66 }67 //68 // Special case for v == 0:69 //70 if (v == 0)71 return -boost::math::detail::cyl_bessel_j_imp<T>(1, x, Tag(), pol);72 //73 // Default case:74 //75 return boost::math::detail::bessel_j_derivative_linear(v, x, Tag(), pol);76}77 78template <class T, class Policy>79inline T sph_bessel_j_prime_imp(unsigned v, T x, const Policy& pol)80{81 static const char* const function = "boost::math::sph_bessel_prime<%1%>(%1%,%1%)";82 //83 // Prevent complex result:84 //85 if (x < 0)86 return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x >= 0.", x, pol);87 //88 // Special case for v == 0:89 //90 if (v == 0)91 return (x == 0) ? boost::math::policies::raise_overflow_error<T>(function, nullptr, pol)92 : static_cast<T>(-boost::math::detail::sph_bessel_j_imp<T>(1, x, pol));93 //94 // Special case for x == 0 and v > 0:95 //96 if (x == 0)97 return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function is indeterminate for this order", x, pol);98 //99 // Default case:100 //101 return boost::math::detail::sph_bessel_j_derivative_linear(v, x, pol);102}103 104template <class T, class Policy>105inline 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_USING109 //110 // Prevent complex result:111 //112 if (x < 0 && floor(v) != v)113 return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x >= 0", x, pol);114 //115 // Special cases for x == 0:116 //117 if (x == 0)118 {119 if (v == 1 || v == -1)120 return static_cast<T>(0.5);121 else if (floor(v) == v || v > 1)122 return 0;123 else 124 return boost::math::policies::raise_domain_error<T>(function, "Got x = %1%, but function is indeterminate for this order", x, pol);125 }126 //127 // Special case for v == 0:128 //129 if (v == 0)130 return boost::math::detail::cyl_bessel_i_imp<T>(1, x, pol);131 //132 // Default case:133 //134 return boost::math::detail::bessel_i_derivative_linear(v, x, pol);135}136 137template <class Tag, class T, class Policy>138inline T cyl_bessel_k_prime_imp(T v, T x, const Policy& pol)139{140 //141 // Prevent complex and indeterminate results:142 //143 if (x <= 0)144 return boost::math::policies::raise_domain_error<T>("boost::math::cyl_bessel_k_prime<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0", x, pol);145 //146 // Special case for v == 0:147 //148 if (v == 0)149 return -boost::math::detail::cyl_bessel_k_imp<T>(1, x, Tag(), pol);150 //151 // Default case:152 //153 return boost::math::detail::bessel_k_derivative_linear(v, x, Tag(), pol);154}155 156template <class Tag, class T, class Policy>157inline T cyl_neumann_prime_imp(T v, T x, const Policy& pol)158{159 BOOST_MATH_STD_USING160 //161 // Prevent complex and indeterminate results:162 //163 if (x <= 0)164 return boost::math::policies::raise_domain_error<T>("boost::math::cyl_neumann_prime<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0", x, pol);165 //166 // Special case for large x: use asymptotic expansion:167 //168 if (boost::math::detail::asymptotic_bessel_derivative_large_x_limit(v, x))169 return boost::math::detail::asymptotic_bessel_y_derivative_large_x_2(v, x, pol);170 //171 // Special case for small x: use Taylor series:172 //173 if (v > 0 && floor(v) != v)174 {175 const T eps = boost::math::policies::get_epsilon<T, Policy>();176 if (log(eps / 2) > v * log((x * x) / (v * 4)))177 return boost::math::detail::bessel_y_derivative_small_z_series(v, x, pol);178 }179 //180 // Special case for v == 0:181 //182 if (v == 0)183 return -boost::math::detail::cyl_neumann_imp<T>(1, x, Tag(), pol);184 //185 // Default case:186 //187 return boost::math::detail::bessel_y_derivative_linear(v, x, Tag(), pol);188}189 190template <class T, class Policy>191inline T sph_neumann_prime_imp(unsigned v, T x, const Policy& pol)192{193 //194 // Prevent complex and indeterminate result:195 //196 if (x <= 0)197 return boost::math::policies::raise_domain_error<T>("boost::math::sph_neumann_prime<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0.", x, pol);198 //199 // Special case for v == 0:200 //201 if (v == 0)202 return -boost::math::detail::sph_neumann_imp<T>(1, x, pol);203 //204 // Default case:205 //206 return boost::math::detail::sph_neumann_derivative_linear(v, x, pol);207}208 209} // namespace detail210 211template <class T1, class T2, class Policy>212inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j_prime(T1 v, T2 x, const Policy& /* pol */)213{214 BOOST_FPU_EXCEPTION_GUARD215 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;216 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;217 typedef typename policies::evaluation<result_type, Policy>::type value_type;218 typedef typename policies::normalise<219 Policy,220 policies::promote_float<false>,221 policies::promote_double<false>,222 policies::discrete_quantile<>,223 policies::assert_undefined<> >::type forwarding_policy;224 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%)");225}226 227template <class T1, class T2>228inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j_prime(T1 v, T2 x)229{230 return cyl_bessel_j_prime(v, x, policies::policy<>());231}232 233template <class T, class Policy>234inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel_prime(unsigned v, T x, const Policy& /* pol */)235{236 BOOST_FPU_EXCEPTION_GUARD237 typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;238 typedef typename policies::evaluation<result_type, Policy>::type value_type;239 typedef typename policies::normalise<240 Policy,241 policies::promote_float<false>,242 policies::promote_double<false>,243 policies::discrete_quantile<>,244 policies::assert_undefined<> >::type forwarding_policy;245 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%)");246}247 248template <class T>249inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel_prime(unsigned v, T x)250{251 return sph_bessel_prime(v, x, policies::policy<>());252}253 254template <class T1, class T2, class Policy>255inline 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_GUARD258 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 269template <class T1, class T2>270inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i_prime(T1 v, T2 x)271{272 return cyl_bessel_i_prime(v, x, policies::policy<>());273}274 275template <class T1, class T2, class Policy>276inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k_prime(T1 v, T2 x, const Policy& /* pol */)277{278 BOOST_FPU_EXCEPTION_GUARD279 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;280 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;281 typedef typename policies::evaluation<result_type, Policy>::type value_type;282 typedef typename policies::normalise<283 Policy,284 policies::promote_float<false>,285 policies::promote_double<false>,286 policies::discrete_quantile<>,287 policies::assert_undefined<> >::type forwarding_policy;288 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%)");289}290 291template <class T1, class T2>292inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k_prime(T1 v, T2 x)293{294 return cyl_bessel_k_prime(v, x, policies::policy<>());295}296 297template <class T1, class T2, class Policy>298inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann_prime(T1 v, T2 x, const Policy& /* pol */)299{300 BOOST_FPU_EXCEPTION_GUARD301 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;302 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;303 typedef typename policies::evaluation<result_type, Policy>::type value_type;304 typedef typename policies::normalise<305 Policy,306 policies::promote_float<false>,307 policies::promote_double<false>,308 policies::discrete_quantile<>,309 policies::assert_undefined<> >::type forwarding_policy;310 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%)");311}312 313template <class T1, class T2>314inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann_prime(T1 v, T2 x)315{316 return cyl_neumann_prime(v, x, policies::policy<>());317}318 319template <class T, class Policy>320inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann_prime(unsigned v, T x, const Policy& /* pol */)321{322 BOOST_FPU_EXCEPTION_GUARD323 typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;324 typedef typename policies::evaluation<result_type, Policy>::type value_type;325 typedef typename policies::normalise<326 Policy,327 policies::promote_float<false>,328 policies::promote_double<false>,329 policies::discrete_quantile<>,330 policies::assert_undefined<> >::type forwarding_policy;331 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%)");332}333 334template <class T>335inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann_prime(unsigned v, T x)336{337 return sph_neumann_prime(v, x, policies::policy<>());338}339 340} // namespace math341} // namespace boost342 343#endif // BOOST_MATH_BESSEL_DERIVATIVES_HPP344