Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-15 08:40:08 

0001 // Copyright John Maddock 2012.
0002 // Copyright Matt Borland 2024.
0003 // Use, modification and distribution are subject to the
0004 // Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt
0006 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0007 
0008 #ifndef BOOST_MATH_AIRY_HPP
0009 #define BOOST_MATH_AIRY_HPP
0010 
0011 #include <boost/math/tools/config.hpp>
0012 #include <boost/math/tools/numeric_limits.hpp>
0013 #include <boost/math/tools/precision.hpp>
0014 #include <boost/math/tools/cstdint.hpp>
0015 #include <boost/math/special_functions/math_fwd.hpp>
0016 #include <boost/math/special_functions/bessel.hpp>
0017 #include <boost/math/special_functions/cbrt.hpp>
0018 #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
0019 #include <boost/math/tools/roots.hpp>
0020 #include <boost/math/policies/error_handling.hpp>
0021 #include <boost/math/constants/constants.hpp>
0022 
0023 namespace boost{ namespace math{
0024 
0025 namespace detail{
0026 
0027 template <class T, class Policy>
0028 BOOST_MATH_GPU_ENABLED T airy_ai_imp(T x, const Policy& pol)
0029 {
0030    BOOST_MATH_STD_USING
0031 
0032    if(x < 0)
0033    {
0034       T p = (-x * sqrt(-x) * 2) / 3;
0035       T v = T(1) / 3;
0036       T j1 = boost::math::cyl_bessel_j(v, p, pol);
0037       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
0038       T ai = sqrt(-x) * (j1 + j2) / 3;
0039       //T bi = sqrt(-x / 3) * (j2 - j1);
0040       return ai;
0041    }
0042    else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
0043    {
0044       T tg = boost::math::tgamma(constants::twothirds<T>(), pol);
0045       T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);
0046       //T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg);
0047       return ai;
0048    }
0049    else
0050    {
0051       T p = 2 * x * sqrt(x) / 3;
0052       T v = T(1) / 3;
0053       //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
0054       //T j2 = boost::math::cyl_bessel_i(v, p, pol);
0055       //
0056       // Note that although we can calculate ai from j1 and j2, the accuracy is horrible
0057       // as we're subtracting two very large values, so use the Bessel K relation instead:
0058       //
0059       T ai = cyl_bessel_k(v, p, pol) * sqrt(x / 3) / boost::math::constants::pi<T>();  //sqrt(x) * (j1 - j2) / 3;
0060       //T bi = sqrt(x / 3) * (j1 + j2);
0061       return ai;
0062    }
0063 }
0064 
0065 template <class T, class Policy>
0066 BOOST_MATH_GPU_ENABLED T airy_bi_imp(T x, const Policy& pol)
0067 {
0068    BOOST_MATH_STD_USING
0069 
0070    if(x < 0)
0071    {
0072       T p = (-x * sqrt(-x) * 2) / 3;
0073       T v = T(1) / 3;
0074       T j1 = boost::math::cyl_bessel_j(v, p, pol);
0075       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
0076       //T ai = sqrt(-x) * (j1 + j2) / 3;
0077       T bi = sqrt(-x / 3) * (j2 - j1);
0078       return bi;
0079    }
0080    else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
0081    {
0082       T tg = boost::math::tgamma(constants::twothirds<T>(), pol);
0083       //T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);
0084       T bi = 1 / (sqrt(boost::math::cbrt(T(3), pol)) * tg);
0085       return bi;
0086    }
0087    else
0088    {
0089       T p = 2 * x * sqrt(x) / 3;
0090       T v = T(1) / 3;
0091       T j1 = boost::math::cyl_bessel_i(-v, p, pol);
0092       T j2 = boost::math::cyl_bessel_i(v, p, pol);
0093       T bi = sqrt(x / 3) * (j1 + j2);
0094       return bi;
0095    }
0096 }
0097 
0098 template <class T, class Policy>
0099 BOOST_MATH_GPU_ENABLED T airy_ai_prime_imp(T x, const Policy& pol)
0100 {
0101    BOOST_MATH_STD_USING
0102 
0103    if(x < 0)
0104    {
0105       T p = (-x * sqrt(-x) * 2) / 3;
0106       T v = T(2) / 3;
0107       T j1 = boost::math::cyl_bessel_j(v, p, pol);
0108       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
0109       T aip = -x * (j1 - j2) / 3;
0110       return aip;
0111    }
0112    else if(fabs(x * x) / 2 < tools::epsilon<T>())
0113    {
0114       T tg = boost::math::tgamma(constants::third<T>(), pol);
0115       T aip = 1 / (boost::math::cbrt(T(3), pol) * tg);
0116       return -aip;
0117    }
0118    else
0119    {
0120       T p = 2 * x * sqrt(x) / 3;
0121       T v = T(2) / 3;
0122       //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
0123       //T j2 = boost::math::cyl_bessel_i(v, p, pol);
0124       //
0125       // Note that although we can calculate ai from j1 and j2, the accuracy is horrible
0126       // as we're subtracting two very large values, so use the Bessel K relation instead:
0127       //
0128       T aip = -cyl_bessel_k(v, p, pol) * x / (boost::math::constants::root_three<T>() * boost::math::constants::pi<T>());
0129       return aip;
0130    }
0131 }
0132 
0133 template <class T, class Policy>
0134 BOOST_MATH_GPU_ENABLED T airy_bi_prime_imp(T x, const Policy& pol)
0135 {
0136    BOOST_MATH_STD_USING
0137 
0138    if(x < 0)
0139    {
0140       T p = (-x * sqrt(-x) * 2) / 3;
0141       T v = T(2) / 3;
0142       T j1 = boost::math::cyl_bessel_j(v, p, pol);
0143       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
0144       T aip = -x * (j1 + j2) / constants::root_three<T>();
0145       return aip;
0146    }
0147    else if(fabs(x * x) / 2 < tools::epsilon<T>())
0148    {
0149       T tg = boost::math::tgamma(constants::third<T>(), pol);
0150       T bip = sqrt(boost::math::cbrt(T(3), pol)) / tg;
0151       return bip;
0152    }
0153    else
0154    {
0155       T p = 2 * x * sqrt(x) / 3;
0156       T v = T(2) / 3;
0157       T j1 = boost::math::cyl_bessel_i(-v, p, pol);
0158       T j2 = boost::math::cyl_bessel_i(v, p, pol);
0159       T aip = x * (j1 + j2) / boost::math::constants::root_three<T>();
0160       return aip;
0161    }
0162 }
0163 
0164 template <class T, class Policy>
0165 BOOST_MATH_GPU_ENABLED T airy_ai_zero_imp(int m, const Policy& pol)
0166 {
0167    BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
0168 
0169    // Handle cases when a negative zero (negative rank) is requested.
0170    if(m < 0)
0171    {
0172       return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%, int)",
0173          "Requested the %1%'th zero, but the rank must be 1 or more !", static_cast<T>(m), pol);
0174    }
0175 
0176    // Handle case when the zero'th zero is requested.
0177    if(m == 0U)
0178    {
0179       return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%,%1%)",
0180         "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);
0181    }
0182 
0183    // Set up the initial guess for the upcoming root-finding.
0184    const T guess_root = boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m, pol);
0185 
0186    // Select the maximum allowed iterations based on the number
0187    // of decimal digits in the numeric type T, being at least 12.
0188    const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F));
0189 
0190    const std::uintmax_t iterations_allowed = static_cast<std::uintmax_t>((std::max)(12, my_digits10 * 2));
0191 
0192    std::uintmax_t iterations_used = iterations_allowed;
0193 
0194    // Use a dynamic tolerance because the roots get closer the higher m gets.
0195    T tolerance;  // LCOV_EXCL_LINE
0196 
0197    if     (m <=   10) { tolerance = T(0.3F); }
0198    else if(m <=  100) { tolerance = T(0.1F); }
0199    else if(m <= 1000) { tolerance = T(0.05F); }
0200    else               { tolerance = T(1) / sqrt(T(m)); }
0201 
0202    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
0203    const T am =
0204       boost::math::tools::newton_raphson_iterate(
0205          boost::math::detail::airy_zero::airy_ai_zero_detail::function_object_ai_and_ai_prime<T, Policy>(pol),
0206          guess_root,
0207          T(guess_root - tolerance),
0208          T(guess_root + tolerance),
0209          policies::digits<T, Policy>(),
0210          iterations_used);
0211 
0212    static_cast<void>(iterations_used);
0213 
0214    return am;
0215 }
0216 
0217 template <class T, class Policy>
0218 BOOST_MATH_GPU_ENABLED T airy_bi_zero_imp(int m, const Policy& pol)
0219 {
0220    BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
0221 
0222    // Handle cases when a negative zero (negative rank) is requested.
0223    if(m < 0)
0224    {
0225       return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%, int)",
0226          "Requested the %1%'th zero, but the rank must 1 or more !", static_cast<T>(m), pol);
0227    }
0228 
0229    // Handle case when the zero'th zero is requested.
0230    if(m == 0U)
0231    {
0232       return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%,%1%)",
0233         "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);
0234    }
0235    // Set up the initial guess for the upcoming root-finding.
0236    const T guess_root = boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m, pol);
0237 
0238    // Select the maximum allowed iterations based on the number
0239    // of decimal digits in the numeric type T, being at least 12.
0240    const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F));
0241 
0242    const std::uintmax_t iterations_allowed = static_cast<std::uintmax_t>((std::max)(12, my_digits10 * 2));
0243 
0244    std::uintmax_t iterations_used = iterations_allowed;
0245 
0246    // Use a dynamic tolerance because the roots get closer the higher m gets.
0247    T tolerance; // LCOV_EXCL_LINE
0248 
0249    if     (m <=   10) { tolerance = T(0.3F); }
0250    else if(m <=  100) { tolerance = T(0.1F); }
0251    else if(m <= 1000) { tolerance = T(0.05F); }
0252    else               { tolerance = T(1) / sqrt(T(m)); }
0253 
0254    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
0255    const T bm =
0256       boost::math::tools::newton_raphson_iterate(
0257          boost::math::detail::airy_zero::airy_bi_zero_detail::function_object_bi_and_bi_prime<T, Policy>(pol),
0258          guess_root,
0259          T(guess_root - tolerance),
0260          T(guess_root + tolerance),
0261          policies::digits<T, Policy>(),
0262          iterations_used);
0263 
0264    static_cast<void>(iterations_used);
0265 
0266    return bm;
0267 }
0268 
0269 } // namespace detail
0270 
0271 template <class T, class Policy>
0272 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai(T x, const Policy&)
0273 {
0274    BOOST_FPU_EXCEPTION_GUARD
0275    typedef typename tools::promote_args<T>::type result_type;
0276    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0277    typedef typename policies::normalise<
0278       Policy, 
0279       policies::promote_float<false>, 
0280       policies::promote_double<false>, 
0281       policies::discrete_quantile<>,
0282       policies::assert_undefined<> >::type forwarding_policy;
0283 
0284    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%)");
0285 }
0286 
0287 template <class T>
0288 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai(T x)
0289 {
0290    return airy_ai(x, policies::policy<>());
0291 }
0292 
0293 template <class T, class Policy>
0294 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi(T x, const Policy&)
0295 {
0296    BOOST_FPU_EXCEPTION_GUARD
0297    typedef typename tools::promote_args<T>::type result_type;
0298    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0299    typedef typename policies::normalise<
0300       Policy, 
0301       policies::promote_float<false>, 
0302       policies::promote_double<false>, 
0303       policies::discrete_quantile<>,
0304       policies::assert_undefined<> >::type forwarding_policy;
0305 
0306    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%)");
0307 }
0308 
0309 template <class T>
0310 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi(T x)
0311 {
0312    return airy_bi(x, policies::policy<>());
0313 }
0314 
0315 template <class T, class Policy>
0316 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai_prime(T x, const Policy&)
0317 {
0318    BOOST_FPU_EXCEPTION_GUARD
0319    typedef typename tools::promote_args<T>::type result_type;
0320    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0321    typedef typename policies::normalise<
0322       Policy, 
0323       policies::promote_float<false>, 
0324       policies::promote_double<false>, 
0325       policies::discrete_quantile<>,
0326       policies::assert_undefined<> >::type forwarding_policy;
0327 
0328    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%)");
0329 }
0330 
0331 template <class T>
0332 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_ai_prime(T x)
0333 {
0334    return airy_ai_prime(x, policies::policy<>());
0335 }
0336 
0337 template <class T, class Policy>
0338 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi_prime(T x, const Policy&)
0339 {
0340    BOOST_FPU_EXCEPTION_GUARD
0341    typedef typename tools::promote_args<T>::type result_type;
0342    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0343    typedef typename policies::normalise<
0344       Policy, 
0345       policies::promote_float<false>, 
0346       policies::promote_double<false>, 
0347       policies::discrete_quantile<>,
0348       policies::assert_undefined<> >::type forwarding_policy;
0349 
0350    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%)");
0351 }
0352 
0353 template <class T>
0354 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type airy_bi_prime(T x)
0355 {
0356    return airy_bi_prime(x, policies::policy<>());
0357 }
0358 
0359 template <class T, class Policy>
0360 BOOST_MATH_GPU_ENABLED inline T airy_ai_zero(int m, const Policy& /*pol*/)
0361 {
0362    BOOST_FPU_EXCEPTION_GUARD
0363    typedef typename policies::evaluation<T, Policy>::type value_type;
0364    typedef typename policies::normalise<
0365       Policy, 
0366       policies::promote_float<false>, 
0367       policies::promote_double<false>, 
0368       policies::discrete_quantile<>,
0369       policies::assert_undefined<> >::type forwarding_policy;
0370 
0371    static_assert(    false == std::numeric_limits<T>::is_specialized
0372                            || (   true  == std::numeric_limits<T>::is_specialized
0373                                && false == std::numeric_limits<T>::is_integer),
0374                            "Airy value type must be a floating-point type.");
0375 
0376    return policies::checked_narrowing_cast<T, Policy>(detail::airy_ai_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_ai_zero<%1%>(unsigned)");
0377 }
0378 
0379 template <class T>
0380 BOOST_MATH_GPU_ENABLED inline T airy_ai_zero(int m)
0381 {
0382    return airy_ai_zero<T>(m, policies::policy<>());
0383 }
0384 
0385 template <class T, class OutputIterator, class Policy>
0386 BOOST_MATH_GPU_ENABLED inline OutputIterator airy_ai_zero(
0387                          int start_index,
0388                          unsigned number_of_zeros,
0389                          OutputIterator out_it,
0390                          const Policy& pol)
0391 {
0392    typedef T result_type;
0393 
0394    static_assert(    false == std::numeric_limits<T>::is_specialized
0395                            || (   true  == std::numeric_limits<T>::is_specialized
0396                                && false == std::numeric_limits<T>::is_integer),
0397                            "Airy value type must be a floating-point type.");
0398 
0399    for(unsigned i = 0; i < number_of_zeros; ++i)
0400    {
0401       *out_it = boost::math::airy_ai_zero<result_type>(start_index + i, pol);
0402       ++out_it;
0403    }
0404    return out_it;
0405 }
0406 
0407 template <class T, class OutputIterator>
0408 BOOST_MATH_GPU_ENABLED inline OutputIterator airy_ai_zero(
0409                          int start_index,
0410                          unsigned number_of_zeros,
0411                          OutputIterator out_it)
0412 {
0413    return airy_ai_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
0414 }
0415 
0416 template <class T, class Policy>
0417 BOOST_MATH_GPU_ENABLED inline T airy_bi_zero(int m, const Policy& /*pol*/)
0418 {
0419    BOOST_FPU_EXCEPTION_GUARD
0420    typedef typename policies::evaluation<T, Policy>::type value_type;
0421    typedef typename policies::normalise<
0422       Policy, 
0423       policies::promote_float<false>, 
0424       policies::promote_double<false>, 
0425       policies::discrete_quantile<>,
0426       policies::assert_undefined<> >::type forwarding_policy;
0427 
0428    static_assert(    false == std::numeric_limits<T>::is_specialized
0429                            || (   true  == std::numeric_limits<T>::is_specialized
0430                                && false == std::numeric_limits<T>::is_integer),
0431                            "Airy value type must be a floating-point type.");
0432 
0433    return policies::checked_narrowing_cast<T, Policy>(detail::airy_bi_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_bi_zero<%1%>(unsigned)");
0434 }
0435 
0436 template <typename T>
0437 BOOST_MATH_GPU_ENABLED inline T airy_bi_zero(int m)
0438 {
0439    return airy_bi_zero<T>(m, policies::policy<>());
0440 }
0441 
0442 template <class T, class OutputIterator, class Policy>
0443 BOOST_MATH_GPU_ENABLED inline OutputIterator airy_bi_zero(
0444                          int start_index,
0445                          unsigned number_of_zeros,
0446                          OutputIterator out_it,
0447                          const Policy& pol)
0448 {
0449    typedef T result_type;
0450 
0451    static_assert(    false == std::numeric_limits<T>::is_specialized
0452                            || (   true  == std::numeric_limits<T>::is_specialized
0453                                && false == std::numeric_limits<T>::is_integer),
0454                            "Airy value type must be a floating-point type.");
0455 
0456    for(unsigned i = 0; i < number_of_zeros; ++i)
0457    {
0458       *out_it = boost::math::airy_bi_zero<result_type>(start_index + i, pol);
0459       ++out_it;
0460    }
0461    return out_it;
0462 }
0463 
0464 template <class T, class OutputIterator>
0465 BOOST_MATH_GPU_ENABLED inline OutputIterator airy_bi_zero(
0466                          int start_index,
0467                          unsigned number_of_zeros,
0468                          OutputIterator out_it)
0469 {
0470    return airy_bi_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
0471 }
0472 
0473 }} // namespaces
0474 
0475 #endif // BOOST_MATH_AIRY_HPP
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!