476 lines · plain
1// Copyright John Maddock 2012.2// Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0.5// (See accompanying file LICENSE_1_0.txt6// or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_AIRY_HPP9#define BOOST_MATH_AIRY_HPP10 11#include <boost/math/tools/config.hpp>12#include <boost/math/tools/numeric_limits.hpp>13#include <boost/math/tools/precision.hpp>14#include <boost/math/tools/cstdint.hpp>15#include <boost/math/special_functions/math_fwd.hpp>16#include <boost/math/special_functions/bessel.hpp>17#include <boost/math/special_functions/cbrt.hpp>18#include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>19#include <boost/math/tools/roots.hpp>20#include <boost/math/policies/error_handling.hpp>21#include <boost/math/constants/constants.hpp>22 23namespace boost{ namespace math{24 25namespace detail{26 27template <class T, class Policy>28BOOST_MATH_GPU_ENABLED T airy_ai_imp(T x, const Policy& pol)29{30 BOOST_MATH_STD_USING31 32 if(x < 0)33 {34 T p = (-x * sqrt(-x) * 2) / 3;35 T v = T(1) / 3;36 T j1 = boost::math::cyl_bessel_j(v, p, pol);37 T j2 = boost::math::cyl_bessel_j(-v, p, pol);38 T ai = sqrt(-x) * (j1 + j2) / 3;39 //T bi = sqrt(-x / 3) * (j2 - j1);40 return ai;41 }42 else if(fabs(x * x * x) / 6 < tools::epsilon<T>())43 {44 T tg = boost::math::tgamma(constants::twothirds<T>(), pol);45 T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);46 //T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg);47 return ai;48 }49 else50 {51 T p = 2 * x * sqrt(x) / 3;52 T v = T(1) / 3;53 //T j1 = boost::math::cyl_bessel_i(-v, p, pol);54 //T j2 = boost::math::cyl_bessel_i(v, p, pol);55 //56 // Note that although we can calculate ai from j1 and j2, the accuracy is horrible57 // as we're subtracting two very large values, so use the Bessel K relation instead:58 //59 T ai = cyl_bessel_k(v, p, pol) * sqrt(x / 3) / boost::math::constants::pi<T>(); //sqrt(x) * (j1 - j2) / 3;60 //T bi = sqrt(x / 3) * (j1 + j2);61 return ai;62 }63}64 65template <class T, class Policy>66BOOST_MATH_GPU_ENABLED T airy_bi_imp(T x, const Policy& pol)67{68 BOOST_MATH_STD_USING69 70 if(x < 0)71 {72 T p = (-x * sqrt(-x) * 2) / 3;73 T v = T(1) / 3;74 T j1 = boost::math::cyl_bessel_j(v, p, pol);75 T j2 = boost::math::cyl_bessel_j(-v, p, pol);76 //T ai = sqrt(-x) * (j1 + j2) / 3;77 T bi = sqrt(-x / 3) * (j2 - j1);78 return bi;79 }80 else if(fabs(x * x * x) / 6 < tools::epsilon<T>())81 {82 T tg = boost::math::tgamma(constants::twothirds<T>(), pol);83 //T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);84 T bi = 1 / (sqrt(boost::math::cbrt(T(3), pol)) * tg);85 return bi;86 }87 else88 {89 T p = 2 * x * sqrt(x) / 3;90 T v = T(1) / 3;91 T j1 = boost::math::cyl_bessel_i(-v, p, pol);92 T j2 = boost::math::cyl_bessel_i(v, p, pol);93 T bi = sqrt(x / 3) * (j1 + j2);94 return bi;95 }96}97 98template <class T, class Policy>99BOOST_MATH_GPU_ENABLED T airy_ai_prime_imp(T x, const Policy& pol)100{101 BOOST_MATH_STD_USING102 103 if(x < 0)104 {105 T p = (-x * sqrt(-x) * 2) / 3;106 T v = T(2) / 3;107 T j1 = boost::math::cyl_bessel_j(v, p, pol);108 T j2 = boost::math::cyl_bessel_j(-v, p, pol);109 T aip = -x * (j1 - j2) / 3;110 return aip;111 }112 else if(fabs(x * x) / 2 < tools::epsilon<T>())113 {114 T tg = boost::math::tgamma(constants::third<T>(), pol);115 T aip = 1 / (boost::math::cbrt(T(3), pol) * tg);116 return -aip;117 }118 else119 {120 T p = 2 * x * sqrt(x) / 3;121 T v = T(2) / 3;122 //T j1 = boost::math::cyl_bessel_i(-v, p, pol);123 //T j2 = boost::math::cyl_bessel_i(v, p, pol);124 //125 // Note that although we can calculate ai from j1 and j2, the accuracy is horrible126 // as we're subtracting two very large values, so use the Bessel K relation instead:127 //128 T aip = -cyl_bessel_k(v, p, pol) * x / (boost::math::constants::root_three<T>() * boost::math::constants::pi<T>());129 return aip;130 }131}132 133template <class T, class Policy>134BOOST_MATH_GPU_ENABLED T airy_bi_prime_imp(T x, const Policy& pol)135{136 BOOST_MATH_STD_USING137 138 if(x < 0)139 {140 T p = (-x * sqrt(-x) * 2) / 3;141 T v = T(2) / 3;142 T j1 = boost::math::cyl_bessel_j(v, p, pol);143 T j2 = boost::math::cyl_bessel_j(-v, p, pol);144 T aip = -x * (j1 + j2) / constants::root_three<T>();145 return aip;146 }147 else if(fabs(x * x) / 2 < tools::epsilon<T>())148 {149 T tg = boost::math::tgamma(constants::third<T>(), pol);150 T bip = sqrt(boost::math::cbrt(T(3), pol)) / tg;151 return bip;152 }153 else154 {155 T p = 2 * x * sqrt(x) / 3;156 T v = T(2) / 3;157 T j1 = boost::math::cyl_bessel_i(-v, p, pol);158 T j2 = boost::math::cyl_bessel_i(v, p, pol);159 T aip = x * (j1 + j2) / boost::math::constants::root_three<T>();160 return aip;161 }162}163 164template <class T, class Policy>165BOOST_MATH_GPU_ENABLED T airy_ai_zero_imp(int m, const Policy& pol)166{167 BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.168 169 // Handle cases when a negative zero (negative rank) is requested.170 if(m < 0)171 {172 return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%, int)",173 "Requested the %1%'th zero, but the rank must be 1 or more !", static_cast<T>(m), pol);174 }175 176 // Handle case when the zero'th zero is requested.177 if(m == 0U)178 {179 return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%,%1%)",180 "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);181 }182 183 // Set up the initial guess for the upcoming root-finding.184 const T guess_root = boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m, pol);185 186 // Select the maximum allowed iterations based on the number187 // of decimal digits in the numeric type T, being at least 12.188 const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F));189 190 const std::uintmax_t iterations_allowed = static_cast<std::uintmax_t>(BOOST_MATH_GPU_SAFE_MAX(12, my_digits10 * 2));191 192 std::uintmax_t iterations_used = iterations_allowed;193 194 // Use a dynamic tolerance because the roots get closer the higher m gets.195 T tolerance; // LCOV_EXCL_LINE196 197 if (m <= 10) { tolerance = T(0.3F); }198 else if(m <= 100) { tolerance = T(0.1F); }199 else if(m <= 1000) { tolerance = T(0.05F); }200 else { tolerance = T(1) / sqrt(T(m)); }201 202 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.203 const T am =204 boost::math::tools::newton_raphson_iterate(205 boost::math::detail::airy_zero::airy_ai_zero_detail::function_object_ai_and_ai_prime<T, Policy>(pol),206 guess_root,207 T(guess_root - tolerance),208 T(guess_root + tolerance),209 policies::digits<T, Policy>(),210 iterations_used);211 212 static_cast<void>(iterations_used);213 214 return am;215}216 217template <class T, class Policy>218BOOST_MATH_GPU_ENABLED T airy_bi_zero_imp(int m, const Policy& pol)219{220 BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.221 222 // Handle cases when a negative zero (negative rank) is requested.223 if(m < 0)224 {225 return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%, int)",226 "Requested the %1%'th zero, but the rank must 1 or more !", static_cast<T>(m), pol);227 }228 229 // Handle case when the zero'th zero is requested.230 if(m == 0U)231 {232 return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%,%1%)",233 "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);234 }235 // Set up the initial guess for the upcoming root-finding.236 const T guess_root = boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m, pol);237 238 // Select the maximum allowed iterations based on the number239 // of decimal digits in the numeric type T, being at least 12.240 const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F));241 242 const std::uintmax_t iterations_allowed = static_cast<std::uintmax_t>(BOOST_MATH_GPU_SAFE_MAX(12, my_digits10 * 2));243 244 std::uintmax_t iterations_used = iterations_allowed;245 246 // Use a dynamic tolerance because the roots get closer the higher m gets.247 T tolerance; // LCOV_EXCL_LINE248 249 if (m <= 10) { tolerance = T(0.3F); }250 else if(m <= 100) { tolerance = T(0.1F); }251 else if(m <= 1000) { tolerance = T(0.05F); }252 else { tolerance = T(1) / sqrt(T(m)); }253 254 // Perform the root-finding using Newton-Raphson iteration from Boost.Math.255 const T bm =256 boost::math::tools::newton_raphson_iterate(257 boost::math::detail::airy_zero::airy_bi_zero_detail::function_object_bi_and_bi_prime<T, Policy>(pol),258 guess_root,259 T(guess_root - tolerance),260 T(guess_root + tolerance),261 policies::digits<T, Policy>(),262 iterations_used);263 264 static_cast<void>(iterations_used);265 266 return bm;267}268 269} // namespace detail270 271template <class T, class Policy>272BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai(T x, const Policy&)273{274 BOOST_FPU_EXCEPTION_GUARD275 typedef typename tools::promote_args<T>::type result_type;276 typedef typename policies::evaluation<result_type, Policy>::type value_type;277 typedef typename policies::normalise<278 Policy, 279 policies::promote_float<false>, 280 policies::promote_double<false>, 281 policies::discrete_quantile<>,282 policies::assert_undefined<> >::type forwarding_policy;283 284 return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_ai_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)");285}286 287template <class T>288BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai(T x)289{290 return airy_ai(x, policies::policy<>());291}292 293template <class T, class Policy>294BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi(T x, const Policy&)295{296 BOOST_FPU_EXCEPTION_GUARD297 typedef typename tools::promote_args<T>::type result_type;298 typedef typename policies::evaluation<result_type, Policy>::type value_type;299 typedef typename policies::normalise<300 Policy, 301 policies::promote_float<false>, 302 policies::promote_double<false>, 303 policies::discrete_quantile<>,304 policies::assert_undefined<> >::type forwarding_policy;305 306 return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_bi_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)");307}308 309template <class T>310BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi(T x)311{312 return airy_bi(x, policies::policy<>());313}314 315template <class T, class Policy>316BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai_prime(T x, const Policy&)317{318 BOOST_FPU_EXCEPTION_GUARD319 typedef typename tools::promote_args<T>::type result_type;320 typedef typename policies::evaluation<result_type, Policy>::type value_type;321 typedef typename policies::normalise<322 Policy, 323 policies::promote_float<false>, 324 policies::promote_double<false>, 325 policies::discrete_quantile<>,326 policies::assert_undefined<> >::type forwarding_policy;327 328 return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_ai_prime_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)");329}330 331template <class T>332BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai_prime(T x)333{334 return airy_ai_prime(x, policies::policy<>());335}336 337template <class T, class Policy>338BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi_prime(T x, const Policy&)339{340 BOOST_FPU_EXCEPTION_GUARD341 typedef typename tools::promote_args<T>::type result_type;342 typedef typename policies::evaluation<result_type, Policy>::type value_type;343 typedef typename policies::normalise<344 Policy, 345 policies::promote_float<false>, 346 policies::promote_double<false>, 347 policies::discrete_quantile<>,348 policies::assert_undefined<> >::type forwarding_policy;349 350 return policies::checked_narrowing_cast<result_type, Policy>(detail::airy_bi_prime_imp<value_type>(static_cast<value_type>(x), forwarding_policy()), "boost::math::airy<%1%>(%1%)");351}352 353template <class T>354BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi_prime(T x)355{356 return airy_bi_prime(x, policies::policy<>());357}358 359template <class T, class Policy>360BOOST_MATH_GPU_ENABLED inline T airy_ai_zero(int m, const Policy& /*pol*/)361{362 BOOST_FPU_EXCEPTION_GUARD363 typedef typename policies::evaluation<T, Policy>::type value_type;364 typedef typename policies::normalise<365 Policy, 366 policies::promote_float<false>, 367 policies::promote_double<false>, 368 policies::discrete_quantile<>,369 policies::assert_undefined<> >::type forwarding_policy;370 371 static_assert( false == std::numeric_limits<T>::is_specialized372 || ( true == std::numeric_limits<T>::is_specialized373 && false == std::numeric_limits<T>::is_integer),374 "Airy value type must be a floating-point type.");375 376 return policies::checked_narrowing_cast<T, Policy>(detail::airy_ai_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_ai_zero<%1%>(unsigned)");377}378 379template <class T>380BOOST_MATH_GPU_ENABLED inline T airy_ai_zero(int m)381{382 return airy_ai_zero<T>(m, policies::policy<>());383}384 385template <class T, class OutputIterator, class Policy>386BOOST_MATH_GPU_ENABLED inline OutputIterator airy_ai_zero(387 int start_index,388 unsigned number_of_zeros,389 OutputIterator out_it,390 const Policy& pol)391{392 typedef T result_type;393 394 static_assert( false == std::numeric_limits<T>::is_specialized395 || ( true == std::numeric_limits<T>::is_specialized396 && false == std::numeric_limits<T>::is_integer),397 "Airy value type must be a floating-point type.");398 399 for(unsigned i = 0; i < number_of_zeros; ++i)400 {401 *out_it = boost::math::airy_ai_zero<result_type>(start_index + i, pol);402 ++out_it;403 }404 return out_it;405}406 407template <class T, class OutputIterator>408BOOST_MATH_GPU_ENABLED inline OutputIterator airy_ai_zero(409 int start_index,410 unsigned number_of_zeros,411 OutputIterator out_it)412{413 return airy_ai_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());414}415 416template <class T, class Policy>417BOOST_MATH_GPU_ENABLED inline T airy_bi_zero(int m, const Policy& /*pol*/)418{419 BOOST_FPU_EXCEPTION_GUARD420 typedef typename policies::evaluation<T, Policy>::type value_type;421 typedef typename policies::normalise<422 Policy, 423 policies::promote_float<false>, 424 policies::promote_double<false>, 425 policies::discrete_quantile<>,426 policies::assert_undefined<> >::type forwarding_policy;427 428 static_assert( false == std::numeric_limits<T>::is_specialized429 || ( true == std::numeric_limits<T>::is_specialized430 && false == std::numeric_limits<T>::is_integer),431 "Airy value type must be a floating-point type.");432 433 return policies::checked_narrowing_cast<T, Policy>(detail::airy_bi_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_bi_zero<%1%>(unsigned)");434}435 436template <typename T>437BOOST_MATH_GPU_ENABLED inline T airy_bi_zero(int m)438{439 return airy_bi_zero<T>(m, policies::policy<>());440}441 442template <class T, class OutputIterator, class Policy>443BOOST_MATH_GPU_ENABLED inline OutputIterator airy_bi_zero(444 int start_index,445 unsigned number_of_zeros,446 OutputIterator out_it,447 const Policy& pol)448{449 typedef T result_type;450 451 static_assert( false == std::numeric_limits<T>::is_specialized452 || ( true == std::numeric_limits<T>::is_specialized453 && false == std::numeric_limits<T>::is_integer),454 "Airy value type must be a floating-point type.");455 456 for(unsigned i = 0; i < number_of_zeros; ++i)457 {458 *out_it = boost::math::airy_bi_zero<result_type>(start_index + i, pol);459 ++out_it;460 }461 return out_it;462}463 464template <class T, class OutputIterator>465BOOST_MATH_GPU_ENABLED inline OutputIterator airy_bi_zero(466 int start_index,467 unsigned number_of_zeros,468 OutputIterator out_it)469{470 return airy_bi_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());471}472 473}} // namespaces474 475#endif // BOOST_MATH_AIRY_HPP476