1 // Copyright (c) 2013 Christopher Kormanyos 2 // Use, modification and distribution are subject to the 3 // Boost Software License, Version 1.0. (See accompanying file 4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 5 // 6 // This work is based on an earlier work: 7 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations", 8 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469 9 // 10 // This header contains implementation details for estimating the zeros 11 // of the Airy functions airy_ai and airy_bi on the negative real axis. 12 // 13 #ifndef BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_ 14 #define BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_ 15 16 #include <boost/math/constants/constants.hpp> 17 #include <boost/math/special_functions/cbrt.hpp> 18 19 namespace boost { namespace math { 20 namespace detail 21 { 22 // Forward declarations of the needed Airy function implementations. 23 template <class T, class Policy> 24 T airy_ai_imp(T x, const Policy& pol); 25 template <class T, class Policy> 26 T airy_bi_imp(T x, const Policy& pol); 27 template <class T, class Policy> 28 T airy_ai_prime_imp(T x, const Policy& pol); 29 template <class T, class Policy> 30 T airy_bi_prime_imp(T x, const Policy& pol); 31 32 namespace airy_zero 33 { 34 template<class T> equation_as_10_4_105(const T & z)35 T equation_as_10_4_105(const T& z) 36 { 37 const T one_over_z (T(1) / z); 38 const T one_over_z_squared(one_over_z * one_over_z); 39 40 const T z_pow_third (boost::math::cbrt(z)); 41 const T z_pow_two_thirds(z_pow_third * z_pow_third); 42 43 // Implement the top line of Eq. 10.4.105. 44 const T fz(z_pow_two_thirds * ((((( + (T(162375596875.0) / 334430208UL) 45 * one_over_z_squared - ( T(108056875.0) / 6967296UL)) 46 * one_over_z_squared + ( T(77125UL) / 82944UL)) 47 * one_over_z_squared - ( T(5U) / 36U)) 48 * one_over_z_squared + ( T(5U) / 48U)) 49 * one_over_z_squared + (1))); 50 51 return fz; 52 } 53 54 namespace airy_ai_zero_detail 55 { 56 template<class T> initial_guess(const int m)57 T initial_guess(const int m) 58 { 59 T guess; 60 61 switch(m) 62 { 63 case 0: { guess = T(0); break; } 64 case 1: { guess = T(-2.33810741045976703849); break; } 65 case 2: { guess = T(-4.08794944413097061664); break; } 66 case 3: { guess = T(-5.52055982809555105913); break; } 67 case 4: { guess = T(-6.78670809007175899878); break; } 68 case 5: { guess = T(-7.94413358712085312314); break; } 69 case 6: { guess = T(-9.02265085334098038016); break; } 70 case 7: { guess = T(-10.0401743415580859306); break; } 71 case 8: { guess = T(-11.0085243037332628932); break; } 72 case 9: { guess = T(-11.9360155632362625170); break; } 73 case 10:{ guess = T(-12.8287767528657572004); break; } 74 default: 75 { 76 const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 1)) / 8); 77 guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t); 78 break; 79 } 80 } 81 82 return guess; 83 } 84 85 template<class T, class Policy> 86 class function_object_ai_and_ai_prime 87 { 88 public: function_object_ai_and_ai_prime(const Policy & pol)89 function_object_ai_and_ai_prime(const Policy& pol) : my_pol(pol) { } 90 operator ()(const T & x) const91 boost::math::tuple<T, T> operator()(const T& x) const 92 { 93 // Return a tuple containing both Ai(x) and Ai'(x). 94 return boost::math::make_tuple( 95 boost::math::detail::airy_ai_imp (x, my_pol), 96 boost::math::detail::airy_ai_prime_imp(x, my_pol)); 97 } 98 99 private: 100 const Policy& my_pol; 101 const function_object_ai_and_ai_prime& operator=(const function_object_ai_and_ai_prime&); 102 }; 103 } // namespace airy_ai_zero_detail 104 105 namespace airy_bi_zero_detail 106 { 107 template<class T> initial_guess(const int m)108 T initial_guess(const int m) 109 { 110 T guess; 111 112 switch(m) 113 { 114 case 0: { guess = T(0); break; } 115 case 1: { guess = T(-1.17371322270912792492); break; } 116 case 2: { guess = T(-3.27109330283635271568); break; } 117 case 3: { guess = T(-4.83073784166201593267); break; } 118 case 4: { guess = T(-6.16985212831025125983); break; } 119 case 5: { guess = T(-7.37676207936776371360); break; } 120 case 6: { guess = T(-8.49194884650938801345); break; } 121 case 7: { guess = T(-9.53819437934623888663); break; } 122 case 8: { guess = T(-10.5299135067053579244); break; } 123 case 9: { guess = T(-11.4769535512787794379); break; } 124 case 10: { guess = T(-12.3864171385827387456); break; } 125 default: 126 { 127 const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 3)) / 8); 128 guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t); 129 break; 130 } 131 } 132 133 return guess; 134 } 135 136 template<class T, class Policy> 137 class function_object_bi_and_bi_prime 138 { 139 public: function_object_bi_and_bi_prime(const Policy & pol)140 function_object_bi_and_bi_prime(const Policy& pol) : my_pol(pol) { } 141 operator ()(const T & x) const142 boost::math::tuple<T, T> operator()(const T& x) const 143 { 144 // Return a tuple containing both Bi(x) and Bi'(x). 145 return boost::math::make_tuple( 146 boost::math::detail::airy_bi_imp (x, my_pol), 147 boost::math::detail::airy_bi_prime_imp(x, my_pol)); 148 } 149 150 private: 151 const Policy& my_pol; 152 const function_object_bi_and_bi_prime& operator=(const function_object_bi_and_bi_prime&); 153 }; 154 } // namespace airy_bi_zero_detail 155 } // namespace airy_zero 156 } // namespace detail 157 } // namespace math 158 } // namespaces boost 159 160 #endif // BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_ 161