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, class Policy>
35 T equation_as_10_4_105(const T& z, const Policy& pol)
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, pol));
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, class Policy>
57 T initial_guess(const int m, const Policy& pol)
58 {
59 T guess;
60
61 switch(m)
62 {
63 case 0:
64 guess = T(0);
65 break;
66 case 1:
67 guess = T(-2.33810741045976703849);
68 break;
69 case 2:
70 guess = T(-4.08794944413097061664);
71 break;
72 case 3:
73 guess = T(-5.52055982809555105913);
74 break;
75 case 4:
76 guess = T(-6.78670809007175899878);
77 break;
78 case 5:
79 guess = T(-7.94413358712085312314);
80 break;
81 case 6:
82 guess = T(-9.02265085334098038016);
83 break;
84 case 7:
85 guess = T(-10.0401743415580859306);
86 break;
87 case 8:
88 guess = T(-11.0085243037332628932);
89 break;
90 case 9:
91 guess = T(-11.9360155632362625170);
92 break;
93 case 10:
94 guess = T(-12.8287767528657572004);
95 break;
96 default:
97 const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 1)) / 8);
98 guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t, pol);
99 break;
100 }
101
102 return guess;
103 }
104
105 template<class T, class Policy>
106 class function_object_ai_and_ai_prime
107 {
108 public:
109 explicit function_object_ai_and_ai_prime(const Policy& pol) : my_pol(pol) { }
110
111 function_object_ai_and_ai_prime(const function_object_ai_and_ai_prime&) = default;
112
113 boost::math::tuple<T, T> operator()(const T& x) const
114 {
115 // Return a tuple containing both Ai(x) and Ai'(x).
116 return boost::math::make_tuple(
117 boost::math::detail::airy_ai_imp (x, my_pol),
118 boost::math::detail::airy_ai_prime_imp(x, my_pol));
119 }
120
121 private:
122 const Policy& my_pol;
123 const function_object_ai_and_ai_prime& operator=(const function_object_ai_and_ai_prime&) = delete;
124 };
125 } // namespace airy_ai_zero_detail
126
127 namespace airy_bi_zero_detail
128 {
129 template<class T, class Policy>
130 T initial_guess(const int m, const Policy& pol)
131 {
132 T guess;
133
134 switch(m)
135 {
136 case 0:
137 guess = T(0);
138 break;
139 case 1:
140 guess = T(-1.17371322270912792492);
141 break;
142 case 2:
143 guess = T(-3.27109330283635271568);
144 break;
145 case 3:
146 guess = T(-4.83073784166201593267);
147 break;
148 case 4:
149 guess = T(-6.16985212831025125983);
150 break;
151 case 5:
152 guess = T(-7.37676207936776371360);
153 break;
154 case 6:
155 guess = T(-8.49194884650938801345);
156 break;
157 case 7:
158 guess = T(-9.53819437934623888663);
159 break;
160 case 8:
161 guess = T(-10.5299135067053579244);
162 break;
163 case 9:
164 guess = T(-11.4769535512787794379);
165 break;
166 case 10:
167 guess = T(-12.3864171385827387456);
168 break;
169 default:
170 const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 3)) / 8);
171 guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t, pol);
172 break;
173 }
174
175 return guess;
176 }
177
178 template<class T, class Policy>
179 class function_object_bi_and_bi_prime
180 {
181 public:
182 explicit function_object_bi_and_bi_prime(const Policy& pol) : my_pol(pol) { }
183
184 function_object_bi_and_bi_prime(const function_object_bi_and_bi_prime&) = default;
185
186 boost::math::tuple<T, T> operator()(const T& x) const
187 {
188 // Return a tuple containing both Bi(x) and Bi'(x).
189 return boost::math::make_tuple(
190 boost::math::detail::airy_bi_imp (x, my_pol),
191 boost::math::detail::airy_bi_prime_imp(x, my_pol));
192 }
193
194 private:
195 const Policy& my_pol;
196 const function_object_bi_and_bi_prime& operator=(const function_object_bi_and_bi_prime&) = delete;
197 };
198 } // namespace airy_bi_zero_detail
199 } // namespace airy_zero
200 } // namespace detail
201 } // namespace math
202 } // namespaces boost
203
204#endif // BOOST_MATH_AIRY_AI_BI_ZERO_2013_01_20_HPP_
205

source code of boost/libs/math/include/boost/math/special_functions/detail/airy_ai_bi_zero.hpp