796 lines · plain
1// Copyright (c) 2007, 2013 John Maddock2// Copyright Christopher Kormanyos 2013.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6//7// This header just defines the function entry points, and adds dispatch8// to the right implementation method. Most of the implementation details9// are in separate headers and copyright Xiaogang Zhang.10//11#ifndef BOOST_MATH_BESSEL_HPP12#define BOOST_MATH_BESSEL_HPP13 14#ifdef _MSC_VER15# pragma once16#endif17 18#include <boost/math/tools/config.hpp>19#include <boost/math/tools/rational.hpp>20#include <boost/math/tools/promotion.hpp>21#include <boost/math/tools/series.hpp>22#include <boost/math/tools/roots.hpp>23#include <boost/math/tools/numeric_limits.hpp>24#include <boost/math/tools/type_traits.hpp>25#include <boost/math/tools/cstdint.hpp>26#include <boost/math/special_functions/detail/bessel_jy.hpp>27#include <boost/math/special_functions/detail/bessel_jn.hpp>28#include <boost/math/special_functions/detail/bessel_yn.hpp>29#include <boost/math/special_functions/detail/bessel_jy_zero.hpp>30#include <boost/math/special_functions/detail/bessel_ik.hpp>31#include <boost/math/special_functions/detail/bessel_i0.hpp>32#include <boost/math/special_functions/detail/bessel_i1.hpp>33#include <boost/math/special_functions/detail/bessel_kn.hpp>34#include <boost/math/special_functions/detail/iconv.hpp>35#include <boost/math/special_functions/sin_pi.hpp>36#include <boost/math/special_functions/cos_pi.hpp>37#include <boost/math/special_functions/sinc.hpp>38#include <boost/math/special_functions/trunc.hpp>39#include <boost/math/special_functions/round.hpp>40#include <boost/math/policies/error_handling.hpp>41#include <boost/math/special_functions/math_fwd.hpp>42 43#ifdef _MSC_VER44# pragma warning(push)45# pragma warning(disable: 6326) // potential comparison of a constant with another constant46#endif47 48namespace boost{ namespace math{49 50namespace detail{51 52template <class T, class Policy>53struct sph_bessel_j_small_z_series_term54{55 typedef T result_type;56 57 BOOST_MATH_GPU_ENABLED sph_bessel_j_small_z_series_term(unsigned v_, T x)58 : N(0), v(v_)59 {60 BOOST_MATH_STD_USING61 mult = x / 2;62 if(v + 3 > max_factorial<T>::value)63 {64 term = v * log(mult) - boost::math::lgamma(v+1+T(0.5f), Policy());65 term = exp(term);66 }67 else68 term = pow(mult, T(v)) / boost::math::tgamma(v+1+T(0.5f), Policy());69 mult *= -mult;70 }71 BOOST_MATH_GPU_ENABLED T operator()()72 {73 T r = term;74 ++N;75 term *= mult / (N * T(N + v + 0.5f));76 return r;77 }78private:79 unsigned N;80 unsigned v;81 T mult;82 T term;83};84 85template <class T, class Policy>86BOOST_MATH_GPU_ENABLED inline T sph_bessel_j_small_z_series(unsigned v, T x, const Policy& pol)87{88 BOOST_MATH_STD_USING // ADL of std names89 sph_bessel_j_small_z_series_term<T, Policy> s(v, x);90 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();91 92 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);93 94 policies::check_series_iterations<T>("boost::math::sph_bessel_j_small_z_series<%1%>(%1%,%1%)", max_iter, pol);95 return result * sqrt(constants::pi<T>() / 4);96}97 98template <class T, class Policy>99BOOST_MATH_GPU_ENABLED T cyl_bessel_j_imp_final(T v, T x, const bessel_no_int_tag&, const Policy& pol)100{101 BOOST_MATH_STD_USING102 103 T result_J, y; // LCOV_EXCL_LINE104 bessel_jy(v, x, &result_J, &y, need_j, pol);105 return result_J;106}107 108// Dispatch funtion to avoid recursion109template <class T, class Policy>110BOOST_MATH_GPU_ENABLED T cyl_bessel_j_imp(T v, T x, const bessel_no_int_tag& t, const Policy& pol)111{112 BOOST_MATH_STD_USING113 114 if(x < 0)115 {116 // better have integer v:117 if (floor(v) == v)118 {119 // LCOV_EXCL_START120 // This branch is hit by multiprecision types only, and is121 // tested by our real_concept tests, but thee are excluded from coverage122 // due to time constraints.123 T r = cyl_bessel_j_imp_final(T(v), T(-x), t, pol);124 if (iround(v, pol) & 1)125 {126 r = -r;127 }128 129 return r;130 // LCOV_EXCL_STOP131 }132 else133 {134 constexpr auto function = "boost::math::bessel_j<%1%>(%1%,%1%)";135 return policies::raise_domain_error<T>(function, "Got x = %1%, but we need x >= 0", x, pol);136 }137 }138 139 return cyl_bessel_j_imp_final(T(v), T(x), t, pol);140}141 142template <class T, class Policy>143BOOST_MATH_GPU_ENABLED inline T cyl_bessel_j_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)144{145 BOOST_MATH_STD_USING // ADL of std names.146 int ival = detail::iconv(v, pol);147 // If v is an integer, use the integer recursion148 // method, both that and Steeds method are O(v):149 if((0 == v - ival))150 {151 return bessel_jn(ival, x, pol);152 }153 return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);154}155 156template <class T, class Policy>157BOOST_MATH_GPU_ENABLED inline T cyl_bessel_j_imp(int v, T x, const bessel_int_tag&, const Policy& pol)158{159 BOOST_MATH_STD_USING160 return bessel_jn(v, x, pol);161}162 163template <class T, class Policy>164BOOST_MATH_GPU_ENABLED inline T sph_bessel_j_imp(unsigned n, T x, const Policy& pol)165{166 BOOST_MATH_STD_USING // ADL of std names167 if(x < 0)168 return policies::raise_domain_error<T>("boost::math::sph_bessel_j<%1%>(%1%,%1%)", "Got x = %1%, but function requires x > 0.", x, pol);169 //170 // Special case, n == 0 resolves down to the sinus cardinal of x:171 //172 if(n == 0)173 return (boost::math::isinf)(x) ? T(0) : boost::math::sinc_pi(x, pol);174 //175 // Special case for x == 0:176 //177 if(x == 0)178 return 0;179 //180 // When x is small we may end up with 0/0, use series evaluation181 // instead, especially as it converges rapidly:182 //183 if(x < 1)184 return sph_bessel_j_small_z_series(n, x, pol);185 //186 // Default case is just a naive evaluation of the definition:187 //188 return sqrt(constants::pi<T>() / (2 * x))189 * cyl_bessel_j_imp(T(T(n)+T(0.5f)), x, bessel_no_int_tag(), pol);190}191 192template <class T, class Policy>193BOOST_MATH_GPU_ENABLED T cyl_bessel_i_imp_final(T v, T x, const Policy& pol)194{195 //196 // This handles all the bessel I functions, note that we don't optimise197 // for integer v, other than the v = 0 or 1 special cases, as Millers198 // algorithm is at least as inefficient as the general case (the general199 // case has better error handling too).200 //201 BOOST_MATH_STD_USING202 constexpr auto function = "boost::math::cyl_bessel_i<%1%>(%1%,%1%)";203 if(x == 0)204 {205 if(v < 0) 206 return floor(v) == v ? static_cast<T>(0) : policies::raise_overflow_error<T>(function, nullptr, pol);207 return (v == 0) ? static_cast<T>(1) : static_cast<T>(0);208 }209 if(v == 0.5f)210 {211 // common special case, note try and avoid overflow in exp(x):212 if(x >= tools::log_max_value<T>())213 {214 T e = exp(x / 2);215 return e * (e / sqrt(2 * x * constants::pi<T>()));216 }217 return sqrt(2 / (x * constants::pi<T>())) * sinh(x);218 }219 if((policies::digits<T, Policy>() <= 113) && (boost::math::numeric_limits<T>::digits <= 113) && (boost::math::numeric_limits<T>::radix == 2))220 {221 if(v == 0)222 {223 return bessel_i0(x);224 }225 if(v == 1)226 {227 return bessel_i1(x);228 }229 }230 if((v > 0) && (x / v < 0.25))231 return bessel_i_small_z_series(v, x, pol);232 T result_I, result_K; // LCOV_EXCL_LINE233 bessel_ik(v, x, &result_I, &result_K, need_i, pol);234 return result_I;235}236 237// Additional dispatch function to get the GPU impls happy238template <class T, class Policy>239BOOST_MATH_GPU_ENABLED T cyl_bessel_i_imp(T v, T x, const Policy& pol)240{241 BOOST_MATH_STD_USING242 constexpr auto function = "boost::math::cyl_bessel_i<%1%>(%1%,%1%)";243 244 if(x < 0)245 {246 // better have integer v:247 if(floor(v) == v)248 {249 T r = cyl_bessel_i_imp_final(T(v), T(-x), pol);250 if(iround(v, pol) & 1)251 {252 r = -r;253 }254 255 return r;256 }257 else258 {259 return policies::raise_domain_error<T>(function, "Got x = %1%, but we need x >= 0", x, pol);260 }261 }262 263 return cyl_bessel_i_imp_final(T(v), T(x), pol);264}265 266template <class T, class Policy>267BOOST_MATH_GPU_ENABLED inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)268{269 constexpr auto function = "boost::math::cyl_bessel_k<%1%>(%1%,%1%)";270 BOOST_MATH_STD_USING271 if(x < 0)272 {273 return policies::raise_domain_error<T>(function, "Got x = %1%, but we need x > 0", x, pol);274 }275 if(x == 0)276 {277 return (v == 0) ? policies::raise_overflow_error<T>(function, nullptr, pol)278 : policies::raise_domain_error<T>(function, "Got x = %1%, but we need x > 0", x, pol);279 }280 T result_I, result_K; // LCOV_EXCL_LINE281 bessel_ik(v, x, &result_I, &result_K, need_k, pol);282 return result_K;283}284 285template <class T, class Policy>286BOOST_MATH_GPU_ENABLED inline T cyl_bessel_k_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)287{288 BOOST_MATH_STD_USING289 if((floor(v) == v))290 {291 return bessel_kn(itrunc(v), x, pol);292 }293 return cyl_bessel_k_imp(v, x, bessel_no_int_tag(), pol);294}295 296template <class T, class Policy>297BOOST_MATH_GPU_ENABLED inline T cyl_bessel_k_imp(int v, T x, const bessel_int_tag&, const Policy& pol)298{299 return bessel_kn(v, x, pol);300}301 302template <class T, class Policy>303BOOST_MATH_GPU_ENABLED inline T cyl_neumann_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)304{305 constexpr auto function = "boost::math::cyl_neumann<%1%>(%1%,%1%)";306 307 BOOST_MATH_INSTRUMENT_VARIABLE(v);308 BOOST_MATH_INSTRUMENT_VARIABLE(x);309 310 if(x <= 0)311 {312 return (v == 0) && (x == 0) ?313 -policies::raise_overflow_error<T>(function, nullptr, pol) // LCOV_EXCL_LINE MP case only here, not tested in code coverage as it takes too long.314 : policies::raise_domain_error<T>(function, "Got x = %1%, but result is complex for x <= 0", x, pol);315 }316 T result_J, y; // LCOV_EXCL_LINE317 bessel_jy(v, x, &result_J, &y, need_y, pol);318 //319 // Post evaluation check for internal overflow during evaluation,320 // can occur when x is small and v is large, in which case the result321 // is -INF:322 //323 if(!(boost::math::isfinite)(y))324 return -policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE defensive programming? Might be dead code now...325 return y;326}327 328template <class T, class Policy>329BOOST_MATH_GPU_ENABLED inline T cyl_neumann_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol)330{331 BOOST_MATH_STD_USING332 333 BOOST_MATH_INSTRUMENT_VARIABLE(v);334 BOOST_MATH_INSTRUMENT_VARIABLE(x);335 336 if(floor(v) == v)337 {338 T r = bessel_yn(itrunc(v, pol), x, pol);339 BOOST_MATH_INSTRUMENT_VARIABLE(r);340 return r;341 }342 T r = cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);343 BOOST_MATH_INSTRUMENT_VARIABLE(r);344 return r;345}346 347template <class T, class Policy>348BOOST_MATH_GPU_ENABLED inline T cyl_neumann_imp(int v, T x, const bessel_int_tag&, const Policy& pol)349{350 return bessel_yn(v, x, pol);351}352 353template <class T, class Policy>354BOOST_MATH_GPU_ENABLED inline T sph_neumann_imp(unsigned v, T x, const Policy& pol)355{356 BOOST_MATH_STD_USING // ADL of std names357 constexpr auto function = "boost::math::sph_neumann<%1%>(%1%,%1%)";358 //359 // Nothing much to do here but check for errors, and360 // evaluate the function's definition directly:361 //362 if(x < 0)363 return policies::raise_domain_error<T>(function, "Got x = %1%, but function requires x > 0.", x, pol);364 365 if(x < 2 * tools::min_value<T>())366 return -policies::raise_overflow_error<T>(function, nullptr, pol);367 368 T result = cyl_neumann_imp(T(T(v)+0.5f), x, bessel_no_int_tag(), pol);369 T tx = sqrt(constants::pi<T>() / (2 * x));370 371 if((tx > 1) && (tools::max_value<T>() / tx < fabs(result)))372 return -policies::raise_overflow_error<T>(function, nullptr, pol);373 374 return result * tx;375}376 377template <class T, class Policy>378BOOST_MATH_GPU_ENABLED inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)379{380 BOOST_MATH_STD_USING // ADL of std names, needed for floor.381 382 constexpr auto function = "boost::math::cyl_bessel_j_zero<%1%>(%1%, int)";383 384 const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);385 386 // Handle non-finite order.387 if (!(boost::math::isfinite)(v) )388 {389 return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);390 }391 392 // Handle negative rank.393 if(m < 0)394 {395 // Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error.396 return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol);397 }398 399 // Get the absolute value of the order.400 const bool order_is_negative = (v < 0);401 const T vv((!order_is_negative) ? v : T(-v));402 403 // Check if the order is very close to zero or very close to an integer.404 const bool order_is_zero = (vv < half_epsilon);405 const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);406 407 if(m == 0)408 {409 if(order_is_zero)410 {411 // The zero'th zero of J0(x) is not defined and requesting it raises a domain error.412 return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of J0, but the rank must be > 0 !", static_cast<T>(m), pol);413 }414 415 // The zero'th zero of Jv(x) for v < 0 is not defined416 // unless the order is a negative integer.417 if(order_is_negative && (!order_is_integer))418 {419 // For non-integer, negative order, requesting the zero'th zero raises a domain error.420 return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Jv for negative, non-integer order, but the rank must be > 0 !", static_cast<T>(m), pol);421 }422 423 // The zero'th zero does exist and its value is zero.424 return T(0);425 }426 427 // Set up the initial guess for the upcoming root-finding.428 // If the order is a negative integer, then use the corresponding429 // positive integer for the order.430 const T guess_root = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess<T, Policy>((order_is_integer ? vv : v), m, pol);431 432 // Select the maximum allowed iterations from the policy.433 boost::math::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();434 435 const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));436 437 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.438 const T jvm =439 boost::math::tools::newton_raphson_iterate(440 boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv_and_jv_prime<T, Policy>((order_is_integer ? vv : v), order_is_zero, pol),441 guess_root,442 T(guess_root - delta_lo),443 T(guess_root + 0.2F),444 policies::digits<T, Policy>(),445 number_of_iterations);446 447 if(number_of_iterations >= policies::get_max_root_iterations<Policy>())448 {449 return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time: Current best guess is %1%", jvm, Policy()); // LCOV_EXCL_LINE450 }451 452 return jvm;453}454 455template <class T, class Policy>456BOOST_MATH_GPU_ENABLED inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)457{458 BOOST_MATH_STD_USING // ADL of std names, needed for floor.459 460 constexpr auto function = "boost::math::cyl_neumann_zero<%1%>(%1%, int)";461 462 // Handle non-finite order.463 if (!(boost::math::isfinite)(v) )464 {465 return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);466 }467 468 // Handle negative rank.469 if(m < 0)470 {471 return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol);472 }473 474 const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);475 476 // Get the absolute value of the order.477 const bool order_is_negative = (v < 0);478 const T vv((!order_is_negative) ? v : T(-v));479 480 const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);481 482 // For negative integers, use reflection to positive integer order.483 if(order_is_negative && order_is_integer)484 return boost::math::detail::cyl_neumann_zero_imp(vv, m, pol);485 486 // Check if the order is very close to a negative half-integer.487 const T delta_half_integer(vv - (floor(vv) + 0.5F));488 489 const bool order_is_negative_half_integer =490 (order_is_negative && ((delta_half_integer > -half_epsilon) && (delta_half_integer < +half_epsilon)));491 492 // The zero'th zero of Yv(x) for v < 0 is not defined493 // unless the order is a negative integer.494 if((m == 0) && (!order_is_negative_half_integer))495 {496 // For non-integer, negative order, requesting the zero'th zero raises a domain error.497 return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Yv for negative, non-half-integer order, but the rank must be > 0 !", static_cast<T>(m), pol);498 }499 500 // For negative half-integers, use the corresponding501 // spherical Bessel function of positive half-integer order.502 if(order_is_negative_half_integer)503 return boost::math::detail::cyl_bessel_j_zero_imp(vv, m, pol);504 505 // Set up the initial guess for the upcoming root-finding.506 // If the order is a negative integer, then use the corresponding507 // positive integer for the order.508 const T guess_root = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess<T, Policy>(v, m, pol);509 510 // Select the maximum allowed iterations from the policy.511 boost::math::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();512 513 const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));514 515 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.516 const T yvm =517 boost::math::tools::newton_raphson_iterate(518 boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol),519 guess_root,520 T(guess_root - delta_lo),521 T(guess_root + 0.2F),522 policies::digits<T, Policy>(),523 number_of_iterations);524 525 if(number_of_iterations >= policies::get_max_root_iterations<Policy>())526 {527 return policies::raise_evaluation_error<T>(function, "Unable to locate root in a reasonable time: Current best guess is %1%", yvm, Policy()); // LCOV_EXCL_LINE528 }529 530 return yvm;531}532 533} // namespace detail534 535template <class T1, class T2, class Policy>536BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_j(T1 v, T2 x, const Policy& /* pol */)537{538 BOOST_FPU_EXCEPTION_GUARD539 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;540 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;541 typedef typename policies::evaluation<result_type, Policy>::type value_type;542 typedef typename policies::normalise<543 Policy,544 policies::promote_float<false>,545 policies::promote_double<false>,546 policies::discrete_quantile<>,547 policies::assert_undefined<> >::type forwarding_policy;548 return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_j<%1%>(%1%,%1%)");549}550 551template <class T1, class T2>552BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_j(T1 v, T2 x)553{554 return cyl_bessel_j(v, x, policies::policy<>());555}556 557template <class T, class Policy>558BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, Policy>::result_type sph_bessel(unsigned v, T x, const Policy& /* pol */)559{560 BOOST_FPU_EXCEPTION_GUARD561 typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;562 typedef typename policies::evaluation<result_type, Policy>::type value_type;563 typedef typename policies::normalise<564 Policy,565 policies::promote_float<false>,566 policies::promote_double<false>,567 policies::discrete_quantile<>,568 policies::assert_undefined<> >::type forwarding_policy;569 return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_bessel_j_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_bessel<%1%>(%1%,%1%)");570}571 572template <class T>573BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_bessel(unsigned v, T x)574{575 return sph_bessel(v, x, policies::policy<>());576}577 578template <class T1, class T2, class Policy>579BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_i(T1 v, T2 x, const Policy& /* pol */)580{581 BOOST_FPU_EXCEPTION_GUARD582 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;583 typedef typename policies::evaluation<result_type, Policy>::type value_type;584 typedef typename policies::normalise<585 Policy,586 policies::promote_float<false>,587 policies::promote_double<false>,588 policies::discrete_quantile<>,589 policies::assert_undefined<> >::type forwarding_policy;590 return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_i_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy()), "boost::math::cyl_bessel_i<%1%>(%1%,%1%)");591}592 593template <class T1, class T2>594BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_i(T1 v, T2 x)595{596 return cyl_bessel_i(v, x, policies::policy<>());597}598 599template <class T1, class T2, class Policy>600BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_bessel_k(T1 v, T2 x, const Policy& /* pol */)601{602 BOOST_FPU_EXCEPTION_GUARD603 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;604 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag128 tag_type;605 typedef typename policies::evaluation<result_type, Policy>::type value_type;606 typedef typename policies::normalise<607 Policy,608 policies::promote_float<false>,609 policies::promote_double<false>,610 policies::discrete_quantile<>,611 policies::assert_undefined<> >::type forwarding_policy;612 return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_k_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_bessel_k<%1%>(%1%,%1%)");613}614 615template <class T1, class T2>616BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_bessel_k(T1 v, T2 x)617{618 return cyl_bessel_k(v, x, policies::policy<>());619}620 621template <class T1, class T2, class Policy>622BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, Policy>::result_type cyl_neumann(T1 v, T2 x, const Policy& /* pol */)623{624 BOOST_FPU_EXCEPTION_GUARD625 typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type;626 typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type;627 typedef typename policies::evaluation<result_type, Policy>::type value_type;628 typedef typename policies::normalise<629 Policy,630 policies::promote_float<false>,631 policies::promote_double<false>,632 policies::discrete_quantile<>,633 policies::assert_undefined<> >::type forwarding_policy;634 return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_imp<value_type>(v, static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::cyl_neumann<%1%>(%1%,%1%)");635}636 637template <class T1, class T2>638BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type cyl_neumann(T1 v, T2 x)639{640 return cyl_neumann(v, x, policies::policy<>());641}642 643template <class T, class Policy>644BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, Policy>::result_type sph_neumann(unsigned v, T x, const Policy& /* pol */)645{646 BOOST_FPU_EXCEPTION_GUARD647 typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;648 typedef typename policies::evaluation<result_type, Policy>::type value_type;649 typedef typename policies::normalise<650 Policy,651 policies::promote_float<false>,652 policies::promote_double<false>,653 policies::discrete_quantile<>,654 policies::assert_undefined<> >::type forwarding_policy;655 return policies::checked_narrowing_cast<result_type, Policy>(detail::sph_neumann_imp<value_type>(v, static_cast<value_type>(x), forwarding_policy()), "boost::math::sph_neumann<%1%>(%1%,%1%)");656}657 658template <class T>659BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann(unsigned v, T x)660{661 return sph_neumann(v, x, policies::policy<>());662}663 664template <class T, class Policy>665BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_bessel_j_zero(T v, int m, const Policy& /* pol */)666{667 BOOST_FPU_EXCEPTION_GUARD668 typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;669 typedef typename policies::evaluation<result_type, Policy>::type value_type;670 typedef typename policies::normalise<671 Policy,672 policies::promote_float<false>,673 policies::promote_double<false>,674 policies::discrete_quantile<>,675 policies::assert_undefined<> >::type forwarding_policy;676 677 static_assert( false == boost::math::numeric_limits<T>::is_specialized678 || ( true == boost::math::numeric_limits<T>::is_specialized679 && false == boost::math::numeric_limits<T>::is_integer),680 "Order must be a floating-point type.");681 682 return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_bessel_j_zero_imp<value_type>(v, m, forwarding_policy()), "boost::math::cyl_bessel_j_zero<%1%>(%1%,%1%)");683}684 685template <class T>686BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_bessel_j_zero(T v, int m)687{688 static_assert( false == boost::math::numeric_limits<T>::is_specialized689 || ( true == boost::math::numeric_limits<T>::is_specialized690 && false == boost::math::numeric_limits<T>::is_integer),691 "Order must be a floating-point type.");692 693 return cyl_bessel_j_zero<T, policies::policy<> >(v, m, policies::policy<>());694}695 696template <class T, class OutputIterator, class Policy>697BOOST_MATH_GPU_ENABLED inline OutputIterator cyl_bessel_j_zero(T v,698 int start_index,699 unsigned number_of_zeros,700 OutputIterator out_it,701 const Policy& pol)702{703 static_assert( false == boost::math::numeric_limits<T>::is_specialized704 || ( true == boost::math::numeric_limits<T>::is_specialized705 && false == boost::math::numeric_limits<T>::is_integer),706 "Order must be a floating-point type.");707 708 for(int i = 0; i < static_cast<int>(number_of_zeros); ++i)709 {710 *out_it = boost::math::cyl_bessel_j_zero(v, start_index + i, pol);711 ++out_it;712 }713 return out_it;714}715 716template <class T, class OutputIterator>717BOOST_MATH_GPU_ENABLED inline OutputIterator cyl_bessel_j_zero(T v,718 int start_index,719 unsigned number_of_zeros,720 OutputIterator out_it)721{722 return cyl_bessel_j_zero(v, start_index, number_of_zeros, out_it, policies::policy<>());723}724 725template <class T, class Policy>726BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, Policy>::result_type cyl_neumann_zero(T v, int m, const Policy& /* pol */)727{728 BOOST_FPU_EXCEPTION_GUARD729 typedef typename detail::bessel_traits<T, T, Policy>::result_type result_type;730 typedef typename policies::evaluation<result_type, Policy>::type value_type;731 typedef typename policies::normalise<732 Policy,733 policies::promote_float<false>,734 policies::promote_double<false>,735 policies::discrete_quantile<>,736 policies::assert_undefined<> >::type forwarding_policy;737 738 static_assert( false == boost::math::numeric_limits<T>::is_specialized739 || ( true == boost::math::numeric_limits<T>::is_specialized740 && false == boost::math::numeric_limits<T>::is_integer),741 "Order must be a floating-point type.");742 743 return policies::checked_narrowing_cast<result_type, Policy>(detail::cyl_neumann_zero_imp<value_type>(v, m, forwarding_policy()), "boost::math::cyl_neumann_zero<%1%>(%1%,%1%)");744}745 746template <class T>747BOOST_MATH_GPU_ENABLED inline typename detail::bessel_traits<T, T, policies::policy<> >::result_type cyl_neumann_zero(T v, int m)748{749 static_assert( false == boost::math::numeric_limits<T>::is_specialized750 || ( true == boost::math::numeric_limits<T>::is_specialized751 && false == boost::math::numeric_limits<T>::is_integer),752 "Order must be a floating-point type.");753 754 return cyl_neumann_zero<T, policies::policy<> >(v, m, policies::policy<>());755}756 757template <class T, class OutputIterator, class Policy>758BOOST_MATH_GPU_ENABLED inline OutputIterator cyl_neumann_zero(T v,759 int start_index,760 unsigned number_of_zeros,761 OutputIterator out_it,762 const Policy& pol)763{764 static_assert( false == boost::math::numeric_limits<T>::is_specialized765 || ( true == boost::math::numeric_limits<T>::is_specialized766 && false == boost::math::numeric_limits<T>::is_integer),767 "Order must be a floating-point type.");768 769 for(int i = 0; i < static_cast<int>(number_of_zeros); ++i)770 {771 *out_it = boost::math::cyl_neumann_zero(v, start_index + i, pol);772 ++out_it;773 }774 return out_it;775}776 777template <class T, class OutputIterator>778BOOST_MATH_GPU_ENABLED inline OutputIterator cyl_neumann_zero(T v,779 int start_index,780 unsigned number_of_zeros,781 OutputIterator out_it)782{783 return cyl_neumann_zero(v, start_index, number_of_zeros, out_it, policies::policy<>());784}785 786} // namespace math787} // namespace boost788 789#ifdef _MSC_VER790# pragma warning(pop)791#endif792 793#endif // BOOST_MATH_BESSEL_HPP794 795 796