]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 _AIRY_AI_BI_ZERO_2013_01_20_HPP_ | |
14 | #define _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> | |
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> | |
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: | |
89 | function_object_ai_and_ai_prime(const Policy pol) : my_pol(pol) { } | |
90 | ||
91 | 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> | |
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: | |
140 | function_object_bi_and_bi_prime(const Policy pol) : my_pol(pol) { } | |
141 | ||
142 | 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 // _AIRY_AI_BI_ZERO_2013_01_20_HPP_ |