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/constants/constants.hpp>
0017   #include <boost/math/special_functions/cbrt.hpp>
0018
0019   namespace boost { namespace math {
0020   namespace detail
0021   {
0022     // Forward declarations of the needed Airy function implementations.
0023     template <class T, class Policy>
0024     T airy_ai_imp(T x, const Policy& pol);
0025     template <class T, class Policy>
0026     T airy_bi_imp(T x, const Policy& pol);
0027     template <class T, class Policy>
0028     T airy_ai_prime_imp(T x, const Policy& pol);
0029     template <class T, class Policy>
0030     T airy_bi_prime_imp(T x, const Policy& pol);
0031
0032     namespace airy_zero
0033     {
0034       template<class T, class Policy>
0035       T equation_as_10_4_105(const T& z, const Policy& pol)
0036       {
0037         const T one_over_z        (T(1) / z);
0038         const T one_over_z_squared(one_over_z * one_over_z);
0039
0040         const T z_pow_third     (boost::math::cbrt(z, pol));
0041         const T z_pow_two_thirds(z_pow_third * z_pow_third);
0042
0043         // Implement the top line of Eq. 10.4.105.
0044         const T fz(z_pow_two_thirds * (((((                     + (T(162375596875.0) / 334430208UL)
0045                                            * one_over_z_squared - (   T(108056875.0) /   6967296UL))
0046                                            * one_over_z_squared + (       T(77125UL) /     82944UL))
0047                                            * one_over_z_squared - (           T(5U)  /        36U))
0048                                            * one_over_z_squared + (           T(5U)  /        48U))
0049                                            * one_over_z_squared + 1));
0050
0051         return fz;
0052       }
0053
0054       namespace airy_ai_zero_detail
0055       {
0056         template<class T, class Policy>
0057         T initial_guess(const int m, const Policy& pol)
0058         {
0059           T guess;
0060
0061           switch(m)
0062           {
0063             case  0:
0064               guess = T(0);
0065               break;
0066             case  1:
0067               guess = T(-2.33810741045976703849);
0068               break;
0069             case  2:
0070               guess = T(-4.08794944413097061664);
0071               break;
0072             case  3:
0073               guess = T(-5.52055982809555105913);
0074               break;
0075             case  4:
0076               guess = T(-6.78670809007175899878);
0077               break;
0078             case  5:
0079               guess = T(-7.94413358712085312314);
0080               break;
0081             case  6:
0082               guess = T(-9.02265085334098038016);
0083               break;
0084             case  7:
0085               guess = T(-10.0401743415580859306);
0086               break;
0087             case  8:
0088               guess = T(-11.0085243037332628932);
0089               break;
0090             case  9:
0091               guess = T(-11.9360155632362625170);
0092               break;
0093             case 10:
0094               guess = T(-12.8287767528657572004);
0095               break;
0096             default:
0097               const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 1)) / 8);
0098               guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t, pol);
0099               break;
0100           }
0101
0102           return guess;
0103         }
0104
0105         template<class T, class Policy>
0106         class function_object_ai_and_ai_prime
0107         {
0108         public:
0109           explicit function_object_ai_and_ai_prime(const Policy& pol) : my_pol(pol) { }
0110
0111           function_object_ai_and_ai_prime(const function_object_ai_and_ai_prime&) = default;
0112
0113           boost::math::tuple<T, T> operator()(const T& x) const
0114           {
0115             // Return a tuple containing both Ai(x) and Ai'(x).
0116             return boost::math::make_tuple(
0117               boost::math::detail::airy_ai_imp      (x, my_pol),
0118               boost::math::detail::airy_ai_prime_imp(x, my_pol));
0119           }
0120
0121         private:
0122           const Policy& my_pol;
0123           const function_object_ai_and_ai_prime& operator=(const function_object_ai_and_ai_prime&) = delete;
0124         };
0125       } // namespace airy_ai_zero_detail
0126
0127       namespace airy_bi_zero_detail
0128       {
0129         template<class T, class Policy>
0130         T initial_guess(const int m, const Policy& pol)
0131         {
0132           T guess;
0133
0134           switch(m)
0135           {
0136             case  0:
0137               guess = T(0);
0138               break;
0139             case  1:
0140               guess = T(-1.17371322270912792492);
0141               break;
0142             case  2:
0143               guess = T(-3.27109330283635271568);
0144               break;
0145             case  3:
0146               guess = T(-4.83073784166201593267);
0147               break;
0148             case  4:
0149               guess = T(-6.16985212831025125983);
0150               break;
0151             case  5:
0152               guess = T(-7.37676207936776371360);
0153               break;
0154             case  6:
0155               guess = T(-8.49194884650938801345);
0156               break;
0157             case  7:
0158               guess = T(-9.53819437934623888663);
0159               break;
0160             case  8:
0161               guess = T(-10.5299135067053579244);
0162               break;
0163             case  9:
0164               guess = T(-11.4769535512787794379);
0165               break;
0166             case 10:
0167               guess = T(-12.3864171385827387456);
0168               break;
0169             default:
0170               const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 3)) / 8);
0171               guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t, pol);
0172               break;
0173           }
0174
0175           return guess;
0176         }
0177
0178         template<class T, class Policy>
0179         class function_object_bi_and_bi_prime
0180         {
0181         public:
0182           explicit function_object_bi_and_bi_prime(const Policy& pol) : my_pol(pol) { }
0183
0184           function_object_bi_and_bi_prime(const function_object_bi_and_bi_prime&) = default;
0185
0186           boost::math::tuple<T, T> operator()(const T& x) const
0187           {
0188             // Return a tuple containing both Bi(x) and Bi'(x).
0189             return boost::math::make_tuple(
0190               boost::math::detail::airy_bi_imp      (x, my_pol),
0191               boost::math::detail::airy_bi_prime_imp(x, my_pol));
0192           }
0193
0194         private:
0195           const Policy& my_pol;
0196           const function_object_bi_and_bi_prime& operator=(const function_object_bi_and_bi_prime&) = delete;
0197         };
0198       } // namespace airy_bi_zero_detail
0199     } // namespace airy_zero
0200   } // namespace detail
0201   } // namespace math
0202   } // namespaces boost
0203
0204 #endif // BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_