Chris@16: // Copyright John Maddock 2012.
Chris@16: // Use, modification and distribution are subject to the
Chris@16: // Boost Software License, Version 1.0.
Chris@16: // (See accompanying file LICENSE_1_0.txt
Chris@16: // or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_AIRY_HPP
Chris@16: #define BOOST_MATH_AIRY_HPP
Chris@16: 
Chris@101: #include <limits>
Chris@101: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/special_functions/bessel.hpp>
Chris@16: #include <boost/math/special_functions/cbrt.hpp>
Chris@16: #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
Chris@16: #include <boost/math/tools/roots.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T airy_ai_imp(T x, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       T p = (-x * sqrt(-x) * 2) / 3;
Chris@16:       T v = T(1) / 3;
Chris@16:       T j1 = boost::math::cyl_bessel_j(v, p, pol);
Chris@16:       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
Chris@16:       T ai = sqrt(-x) * (j1 + j2) / 3;
Chris@16:       //T bi = sqrt(-x / 3) * (j2 - j1);
Chris@16:       return ai;
Chris@16:    }
Chris@16:    else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
Chris@16:    {
Chris@16:       T tg = boost::math::tgamma(constants::twothirds<T>(), pol);
Chris@16:       T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);
Chris@16:       //T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg);
Chris@16:       return ai;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       T p = 2 * x * sqrt(x) / 3;
Chris@16:       T v = T(1) / 3;
Chris@16:       //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
Chris@16:       //T j2 = boost::math::cyl_bessel_i(v, p, pol);
Chris@16:       //
Chris@16:       // Note that although we can calculate ai from j1 and j2, the accuracy is horrible
Chris@16:       // as we're subtracting two very large values, so use the Bessel K relation instead:
Chris@16:       //
Chris@16:       T ai = cyl_bessel_k(v, p, pol) * sqrt(x / 3) / boost::math::constants::pi<T>();  //sqrt(x) * (j1 - j2) / 3;
Chris@16:       //T bi = sqrt(x / 3) * (j1 + j2);
Chris@16:       return ai;
Chris@16:    }
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T airy_bi_imp(T x, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       T p = (-x * sqrt(-x) * 2) / 3;
Chris@16:       T v = T(1) / 3;
Chris@16:       T j1 = boost::math::cyl_bessel_j(v, p, pol);
Chris@16:       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
Chris@16:       //T ai = sqrt(-x) * (j1 + j2) / 3;
Chris@16:       T bi = sqrt(-x / 3) * (j2 - j1);
Chris@16:       return bi;
Chris@16:    }
Chris@16:    else if(fabs(x * x * x) / 6 < tools::epsilon<T>())
Chris@16:    {
Chris@16:       T tg = boost::math::tgamma(constants::twothirds<T>(), pol);
Chris@16:       //T ai = 1 / (pow(T(3), constants::twothirds<T>()) * tg);
Chris@16:       T bi = 1 / (sqrt(boost::math::cbrt(T(3))) * tg);
Chris@16:       return bi;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       T p = 2 * x * sqrt(x) / 3;
Chris@16:       T v = T(1) / 3;
Chris@16:       T j1 = boost::math::cyl_bessel_i(-v, p, pol);
Chris@16:       T j2 = boost::math::cyl_bessel_i(v, p, pol);
Chris@16:       T bi = sqrt(x / 3) * (j1 + j2);
Chris@16:       return bi;
Chris@16:    }
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T airy_ai_prime_imp(T x, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       T p = (-x * sqrt(-x) * 2) / 3;
Chris@16:       T v = T(2) / 3;
Chris@16:       T j1 = boost::math::cyl_bessel_j(v, p, pol);
Chris@16:       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
Chris@16:       T aip = -x * (j1 - j2) / 3;
Chris@16:       return aip;
Chris@16:    }
Chris@16:    else if(fabs(x * x) / 2 < tools::epsilon<T>())
Chris@16:    {
Chris@16:       T tg = boost::math::tgamma(constants::third<T>(), pol);
Chris@16:       T aip = 1 / (boost::math::cbrt(T(3)) * tg);
Chris@16:       return -aip;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       T p = 2 * x * sqrt(x) / 3;
Chris@16:       T v = T(2) / 3;
Chris@16:       //T j1 = boost::math::cyl_bessel_i(-v, p, pol);
Chris@16:       //T j2 = boost::math::cyl_bessel_i(v, p, pol);
Chris@16:       //
Chris@16:       // Note that although we can calculate ai from j1 and j2, the accuracy is horrible
Chris@16:       // as we're subtracting two very large values, so use the Bessel K relation instead:
Chris@16:       //
Chris@16:       T aip = -cyl_bessel_k(v, p, pol) * x / (boost::math::constants::root_three<T>() * boost::math::constants::pi<T>());
Chris@16:       return aip;
Chris@16:    }
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T airy_bi_prime_imp(T x, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       T p = (-x * sqrt(-x) * 2) / 3;
Chris@16:       T v = T(2) / 3;
Chris@16:       T j1 = boost::math::cyl_bessel_j(v, p, pol);
Chris@16:       T j2 = boost::math::cyl_bessel_j(-v, p, pol);
Chris@16:       T aip = -x * (j1 + j2) / constants::root_three<T>();
Chris@16:       return aip;
Chris@16:    }
Chris@16:    else if(fabs(x * x) / 2 < tools::epsilon<T>())
Chris@16:    {
Chris@16:       T tg = boost::math::tgamma(constants::third<T>(), pol);
Chris@16:       T bip = sqrt(boost::math::cbrt(T(3))) / tg;
Chris@16:       return bip;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       T p = 2 * x * sqrt(x) / 3;
Chris@16:       T v = T(2) / 3;
Chris@16:       T j1 = boost::math::cyl_bessel_i(-v, p, pol);
Chris@16:       T j2 = boost::math::cyl_bessel_i(v, p, pol);
Chris@16:       T aip = x * (j1 + j2) / boost::math::constants::root_three<T>();
Chris@16:       return aip;
Chris@16:    }
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@101: T airy_ai_zero_imp(int m, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
Chris@16: 
Chris@101:    // Handle cases when a negative zero (negative rank) is requested.
Chris@101:    if(m < 0)
Chris@101:    {
Chris@101:       return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%, int)",
Chris@101:          "Requested the %1%'th zero, but the rank must be 1 or more !", static_cast<T>(m), pol);
Chris@101:    }
Chris@101: 
Chris@16:    // Handle case when the zero'th zero is requested.
Chris@16:    if(m == 0U)
Chris@16:    {
Chris@16:       return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%,%1%)",
Chris@16:         "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);
Chris@16:    }
Chris@16: 
Chris@16:    // Set up the initial guess for the upcoming root-finding.
Chris@16:    const T guess_root = boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m);
Chris@16: 
Chris@101:    // Select the maximum allowed iterations based on the number
Chris@101:    // of decimal digits in the numeric type T, being at least 12.
Chris@101:    const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F));
Chris@16: 
Chris@101:    const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2));
Chris@101: 
Chris@101:    boost::uintmax_t iterations_used = iterations_allowed;
Chris@16: 
Chris@16:    // Use a dynamic tolerance because the roots get closer the higher m gets.
Chris@16:    T tolerance;
Chris@16: 
Chris@101:    if     (m <=   10) { tolerance = T(0.3F); }
Chris@101:    else if(m <=  100) { tolerance = T(0.1F); }
Chris@101:    else if(m <= 1000) { tolerance = T(0.05F); }
Chris@101:    else               { tolerance = T(1) / sqrt(T(m)); }
Chris@16: 
Chris@16:    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
Chris@16:    const T am =
Chris@16:       boost::math::tools::newton_raphson_iterate(
Chris@16:          boost::math::detail::airy_zero::airy_ai_zero_detail::function_object_ai_and_ai_prime<T, Policy>(pol),
Chris@16:          guess_root,
Chris@16:          T(guess_root - tolerance),
Chris@16:          T(guess_root + tolerance),
Chris@101:          policies::digits<T, Policy>(),
Chris@101:          iterations_used);
Chris@16: 
Chris@101:    static_cast<void>(iterations_used);
Chris@16: 
Chris@16:    return am;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@101: T airy_bi_zero_imp(int m, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
Chris@16: 
Chris@101:    // Handle cases when a negative zero (negative rank) is requested.
Chris@101:    if(m < 0)
Chris@101:    {
Chris@101:       return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%, int)",
Chris@101:          "Requested the %1%'th zero, but the rank must 1 or more !", static_cast<T>(m), pol);
Chris@101:    }
Chris@101: 
Chris@16:    // Handle case when the zero'th zero is requested.
Chris@16:    if(m == 0U)
Chris@16:    {
Chris@16:       return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%,%1%)",
Chris@16:         "The requested rank of the zero is %1%, but must be 1 or more !", static_cast<T>(m), pol);
Chris@16:    }
Chris@16:    // Set up the initial guess for the upcoming root-finding.
Chris@16:    const T guess_root = boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m);
Chris@16: 
Chris@101:    // Select the maximum allowed iterations based on the number
Chris@101:    // of decimal digits in the numeric type T, being at least 12.
Chris@101:    const int my_digits10 = static_cast<int>(static_cast<float>(policies::digits<T, Policy>() * 0.301F));
Chris@16: 
Chris@101:    const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2));
Chris@101: 
Chris@101:    boost::uintmax_t iterations_used = iterations_allowed;
Chris@16: 
Chris@16:    // Use a dynamic tolerance because the roots get closer the higher m gets.
Chris@16:    T tolerance;
Chris@16: 
Chris@101:    if     (m <=   10) { tolerance = T(0.3F); }
Chris@101:    else if(m <=  100) { tolerance = T(0.1F); }
Chris@101:    else if(m <= 1000) { tolerance = T(0.05F); }
Chris@101:    else               { tolerance = T(1) / sqrt(T(m)); }
Chris@16: 
Chris@16:    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
Chris@16:    const T bm =
Chris@16:       boost::math::tools::newton_raphson_iterate(
Chris@16:          boost::math::detail::airy_zero::airy_bi_zero_detail::function_object_bi_and_bi_prime<T, Policy>(pol),
Chris@16:          guess_root,
Chris@16:          T(guess_root - tolerance),
Chris@16:          T(guess_root + tolerance),
Chris@101:          policies::digits<T, Policy>(),
Chris@101:          iterations_used);
Chris@16: 
Chris@101:    static_cast<void>(iterations_used);
Chris@16: 
Chris@16:    return bm;
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type airy_ai(T x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type airy_ai(T x)
Chris@16: {
Chris@16:    return airy_ai(x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type airy_bi(T x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type airy_bi(T x)
Chris@16: {
Chris@16:    return airy_bi(x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type airy_ai_prime(T x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type airy_ai_prime(T x)
Chris@16: {
Chris@16:    return airy_ai_prime(x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type airy_bi_prime(T x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type airy_bi_prime(T x)
Chris@16: {
Chris@16:    return airy_bi_prime(x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@101: inline T airy_ai_zero(int m, const Policy& /*pol*/)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename policies::evaluation<T, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@101: 
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Airy value type must be a floating-point type.");
Chris@101: 
Chris@16:    return policies::checked_narrowing_cast<T, Policy>(detail::airy_ai_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_ai_zero<%1%>(unsigned)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@101: inline T airy_ai_zero(int m)
Chris@16: {
Chris@16:    return airy_ai_zero<T>(m, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator, class Policy>
Chris@16: inline OutputIterator airy_ai_zero(
Chris@101:                          int start_index,
Chris@16:                          unsigned number_of_zeros,
Chris@16:                          OutputIterator out_it,
Chris@16:                          const Policy& pol)
Chris@16: {
Chris@16:    typedef T result_type;
Chris@101: 
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Airy value type must be a floating-point type.");
Chris@16: 
Chris@16:    for(unsigned i = 0; i < number_of_zeros; ++i)
Chris@16:    {
Chris@16:       *out_it = boost::math::airy_ai_zero<result_type>(start_index + i, pol);
Chris@16:       ++out_it;
Chris@16:    }
Chris@16:    return out_it;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator>
Chris@16: inline OutputIterator airy_ai_zero(
Chris@101:                          int start_index,
Chris@16:                          unsigned number_of_zeros,
Chris@16:                          OutputIterator out_it)
Chris@16: {
Chris@16:    return airy_ai_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@101: inline T airy_bi_zero(int m, const Policy& /*pol*/)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename policies::evaluation<T, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@101: 
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Airy value type must be a floating-point type.");
Chris@101: 
Chris@16:    return policies::checked_narrowing_cast<T, Policy>(detail::airy_bi_zero_imp<value_type>(m, forwarding_policy()), "boost::math::airy_bi_zero<%1%>(unsigned)");
Chris@16: }
Chris@16: 
Chris@16: template <typename T>
Chris@101: inline T airy_bi_zero(int m)
Chris@16: {
Chris@16:    return airy_bi_zero<T>(m, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator, class Policy>
Chris@16: inline OutputIterator airy_bi_zero(
Chris@101:                          int start_index,
Chris@16:                          unsigned number_of_zeros,
Chris@16:                          OutputIterator out_it,
Chris@16:                          const Policy& pol)
Chris@16: {
Chris@16:    typedef T result_type;
Chris@101: 
Chris@101:    BOOST_STATIC_ASSERT_MSG(    false == std::numeric_limits<T>::is_specialized
Chris@101:                            || (   true  == std::numeric_limits<T>::is_specialized
Chris@101:                                && false == std::numeric_limits<T>::is_integer),
Chris@101:                            "Airy value type must be a floating-point type.");
Chris@16: 
Chris@16:    for(unsigned i = 0; i < number_of_zeros; ++i)
Chris@16:    {
Chris@16:       *out_it = boost::math::airy_bi_zero<result_type>(start_index + i, pol);
Chris@16:       ++out_it;
Chris@16:    }
Chris@16:    return out_it;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class OutputIterator>
Chris@16: inline OutputIterator airy_bi_zero(
Chris@101:                          int start_index,
Chris@16:                          unsigned number_of_zeros,
Chris@16:                          OutputIterator out_it)
Chris@16: {
Chris@16:    return airy_bi_zero<T>(start_index, number_of_zeros, out_it, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: }} // namespaces
Chris@16: 
Chris@16: #endif // BOOST_MATH_AIRY_HPP