special_functions/detail/airy_ai_bi_zero.hpp

0001 //  Copyright (c) 2013 Christopher Kormanyos
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005 //
0006 // This work is based on an earlier work:
0007 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
0008 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
0009 //
0010 // This header contains implementation details for estimating the zeros
0011 // of the Airy functions airy_ai and airy_bi on the negative real axis.
0012 //
0013 #ifndef BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_
0014   #define BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_
0015
0016   #include <boost/math/tools/config.hpp>
0017   #include <boost/math/tools/tuple.hpp>
0018   #include <boost/math/constants/constants.hpp>
0019   #include <boost/math/special_functions/cbrt.hpp>
0020
0021   namespace boost { namespace math {
0022   namespace detail
0023   {
0024     // Forward declarations of the needed Airy function implementations.
0025     template <class T, class Policy>
0026     BOOST_MATH_GPU_ENABLED T airy_ai_imp(T x, const Policy& pol);
0027     template <class T, class Policy>
0028     BOOST_MATH_GPU_ENABLED T airy_bi_imp(T x, const Policy& pol);
0029     template <class T, class Policy>
0030     BOOST_MATH_GPU_ENABLED T airy_ai_prime_imp(T x, const Policy& pol);
0031     template <class T, class Policy>
0032     BOOST_MATH_GPU_ENABLED T airy_bi_prime_imp(T x, const Policy& pol);
0033
0034     namespace airy_zero
0035     {
0036       template<class T, class Policy>
0037       BOOST_MATH_GPU_ENABLED T equation_as_10_4_105(const T& z, const Policy& pol)
0038       {
0039         const T one_over_z        (T(1) / z);
0040         const T one_over_z_squared(one_over_z * one_over_z);
0041
0042         const T z_pow_third     (boost::math::cbrt(z, pol));
0043         const T z_pow_two_thirds(z_pow_third * z_pow_third);
0044
0045         // Implement the top line of Eq. 10.4.105.
0046         const T fz(z_pow_two_thirds * (((((                     + (T(162375596875.0) / 334430208UL)
0047                                            * one_over_z_squared - (   T(108056875.0) /   6967296UL))
0048                                            * one_over_z_squared + (       T(77125UL) /     82944UL))
0049                                            * one_over_z_squared - (           T(5U)  /        36U))
0050                                            * one_over_z_squared + (           T(5U)  /        48U))
0051                                            * one_over_z_squared + 1));
0052
0053         return fz;
0054       }
0055
0056       namespace airy_ai_zero_detail
0057       {
0058         template<class T, class Policy>
0059         BOOST_MATH_GPU_ENABLED T initial_guess(const int m, const Policy& pol)
0060         {
0061           T guess;
0062
0063           switch(m)
0064           {
0065             case  0:
0066               guess = T(0);
0067               break;
0068             case  1:
0069               guess = T(-2.33810741045976703849);
0070               break;
0071             case  2:
0072               guess = T(-4.08794944413097061664);
0073               break;
0074             case  3:
0075               guess = T(-5.52055982809555105913);
0076               break;
0077             case  4:
0078               guess = T(-6.78670809007175899878);
0079               break;
0080             case  5:
0081               guess = T(-7.94413358712085312314);
0082               break;
0083             case  6:
0084               guess = T(-9.02265085334098038016);
0085               break;
0086             case  7:
0087               guess = T(-10.0401743415580859306);
0088               break;
0089             case  8:
0090               guess = T(-11.0085243037332628932);
0091               break;
0092             case  9:
0093               guess = T(-11.9360155632362625170);
0094               break;
0095             case 10:
0096               guess = T(-12.8287767528657572004);
0097               break;
0098             default:
0099               const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 1)) / 8);
0100               guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t, pol);
0101               break;
0102           }
0103
0104           return guess;
0105         }
0106
0107         template<class T, class Policy>
0108         class function_object_ai_and_ai_prime
0109         {
0110         public:
0111           BOOST_MATH_GPU_ENABLED explicit function_object_ai_and_ai_prime(const Policy& pol) : my_pol(pol) { }
0112
0113           #ifdef BOOST_MATH_ENABLE_CUDA
0114           #  pragma nv_diag_suppress 20012
0115           #endif
0116
0117           BOOST_MATH_GPU_ENABLED function_object_ai_and_ai_prime(const function_object_ai_and_ai_prime&) = default;
0118
0119           #ifdef BOOST_MATH_ENABLE_CUDA
0120           #  pragma nv_diag_default 20012
0121           #endif
0122
0123           BOOST_MATH_GPU_ENABLED boost::math::tuple<T, T> operator()(const T& x) const
0124           {
0125             // Return a tuple containing both Ai(x) and Ai'(x).
0126             return boost::math::make_tuple(
0127               boost::math::detail::airy_ai_imp      (x, my_pol),
0128               boost::math::detail::airy_ai_prime_imp(x, my_pol));
0129           }
0130
0131         private:
0132           const Policy& my_pol;
0133           const function_object_ai_and_ai_prime& operator=(const function_object_ai_and_ai_prime&) = delete;
0134         };
0135       } // namespace airy_ai_zero_detail
0136
0137       namespace airy_bi_zero_detail
0138       {
0139         template<class T, class Policy>
0140         BOOST_MATH_GPU_ENABLED T initial_guess(const int m, const Policy& pol)
0141         {
0142           T guess;
0143
0144           switch(m)
0145           {
0146             case  0:
0147               guess = T(0);
0148               break;
0149             case  1:
0150               guess = T(-1.17371322270912792492);
0151               break;
0152             case  2:
0153               guess = T(-3.27109330283635271568);
0154               break;
0155             case  3:
0156               guess = T(-4.83073784166201593267);
0157               break;
0158             case  4:
0159               guess = T(-6.16985212831025125983);
0160               break;
0161             case  5:
0162               guess = T(-7.37676207936776371360);
0163               break;
0164             case  6:
0165               guess = T(-8.49194884650938801345);
0166               break;
0167             case  7:
0168               guess = T(-9.53819437934623888663);
0169               break;
0170             case  8:
0171               guess = T(-10.5299135067053579244);
0172               break;
0173             case  9:
0174               guess = T(-11.4769535512787794379);
0175               break;
0176             case 10:
0177               guess = T(-12.3864171385827387456);
0178               break;
0179             default:
0180               const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 3)) / 8);
0181               guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t, pol);
0182               break;
0183           }
0184
0185           return guess;
0186         }
0187
0188         template<class T, class Policy>
0189         class function_object_bi_and_bi_prime
0190         {
0191         public:
0192           BOOST_MATH_GPU_ENABLED explicit function_object_bi_and_bi_prime(const Policy& pol) : my_pol(pol) { }
0193
0194           #ifdef BOOST_MATH_ENABLE_CUDA
0195           #  pragma nv_diag_suppress 20012
0196           #endif
0197
0198           BOOST_MATH_GPU_ENABLED function_object_bi_and_bi_prime(const function_object_bi_and_bi_prime&) = default;
0199
0200           #ifdef BOOST_MATH_ENABLE_CUDA
0201           #  pragma nv_diag_default 20012
0202           #endif
0203
0204           BOOST_MATH_GPU_ENABLED boost::math::tuple<T, T> operator()(const T& x) const
0205           {
0206             // Return a tuple containing both Bi(x) and Bi'(x).
0207             return boost::math::make_tuple(
0208               boost::math::detail::airy_bi_imp      (x, my_pol),
0209               boost::math::detail::airy_bi_prime_imp(x, my_pol));
0210           }
0211
0212         private:
0213           const Policy& my_pol;
0214           const function_object_bi_and_bi_prime& operator=(const function_object_bi_and_bi_prime&) = delete;
0215         };
0216       } // namespace airy_bi_zero_detail
0217     } // namespace airy_zero
0218   } // namespace detail
0219   } // namespace math
0220   } // namespaces boost
0221
0222 #endif // BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_