math/special_functions/airy.hpp

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