]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | // (C) Copyright Anton Bikineev 2014 |
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 | #ifndef BOOST_MATH_TOOLS_RECURRENCE_HPP_ | |
7 | #define BOOST_MATH_TOOLS_RECURRENCE_HPP_ | |
8 | ||
9 | #include <boost/math/tools/config.hpp> | |
10 | #include <boost/math/tools/precision.hpp> | |
11 | #include <boost/math/tools/tuple.hpp> | |
12 | #include <boost/math/tools/fraction.hpp> | |
13 | ||
14 | #ifdef BOOST_NO_CXX11_HDR_TUPLE | |
15 | #error "This header requires C++11 support" | |
16 | #endif | |
17 | ||
18 | ||
19 | namespace boost { | |
20 | namespace math { | |
21 | namespace tools { | |
22 | namespace detail{ | |
23 | ||
24 | // | |
25 | // Function ratios directly from recurrence relations: | |
26 | // H. Shintan, Note on Miller's recurrence algorithm, J. Sci. Hiroshima Univ. Ser. A-I | |
27 | // Math., 29 (1965), pp. 121 - 133. | |
28 | // and: | |
29 | // COMPUTATIONAL ASPECTS OF THREE-TERM RECURRENCE RELATIONS | |
30 | // WALTER GAUTSCHI | |
31 | // SIAM REVIEW Vol. 9, No. 1, January, 1967 | |
32 | // | |
33 | template <class Recurrence> | |
34 | struct function_ratio_from_backwards_recurrence_fraction | |
35 | { | |
36 | typedef typename boost::remove_reference<decltype(boost::math::get<0>(std::declval<Recurrence&>()(0)))>::type value_type; | |
37 | typedef std::pair<value_type, value_type> result_type; | |
38 | function_ratio_from_backwards_recurrence_fraction(const Recurrence& r) : r(r), k(0) {} | |
39 | ||
40 | result_type operator()() | |
41 | { | |
42 | value_type a, b, c; | |
43 | boost::math::tie(a, b, c) = r(k); | |
44 | ++k; | |
45 | // an and bn defined as per Gauchi 1.16, not the same | |
46 | // as the usual continued fraction a' and b's. | |
47 | value_type bn = a / c; | |
48 | value_type an = b / c; | |
49 | return result_type(-bn, an); | |
50 | } | |
51 | ||
52 | private: | |
53 | function_ratio_from_backwards_recurrence_fraction operator=(const function_ratio_from_backwards_recurrence_fraction&); | |
54 | ||
55 | Recurrence r; | |
56 | int k; | |
57 | }; | |
58 | ||
59 | template <class R, class T> | |
60 | struct recurrence_reverser | |
61 | { | |
62 | recurrence_reverser(const R& r) : r(r) {} | |
63 | boost::math::tuple<T, T, T> operator()(int i) | |
64 | { | |
65 | using std::swap; | |
66 | boost::math::tuple<T, T, T> t = r(-i); | |
67 | swap(boost::math::get<0>(t), boost::math::get<2>(t)); | |
68 | return t; | |
69 | } | |
70 | R r; | |
71 | }; | |
72 | ||
73 | template <class Recurrence> | |
74 | struct recurrence_offsetter | |
75 | { | |
76 | typedef decltype(std::declval<Recurrence&>()(0)) result_type; | |
77 | recurrence_offsetter(Recurrence const& rr, int offset) : r(rr), k(offset) {} | |
78 | result_type operator()(int i) | |
79 | { | |
80 | return r(i + k); | |
81 | } | |
82 | private: | |
83 | Recurrence r; | |
84 | int k; | |
85 | }; | |
86 | ||
87 | ||
88 | ||
89 | } // namespace detail | |
90 | ||
91 | // | |
92 | // Given a stable backwards recurrence relation: | |
93 | // a f_n-1 + b f_n + c f_n+1 = 0 | |
94 | // returns the ratio f_n / f_n-1 | |
95 | // | |
96 | // Recurrence: a functor that returns a tuple of the factors (a,b,c). | |
97 | // factor: Convergence criteria, should be no less than machine epsilon. | |
98 | // max_iter: Maximum iterations to use solving the continued fraction. | |
99 | // | |
100 | template <class Recurrence, class T> | |
101 | T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter) | |
102 | { | |
103 | detail::function_ratio_from_backwards_recurrence_fraction<Recurrence> f(r); | |
104 | return boost::math::tools::continued_fraction_a(f, factor, max_iter); | |
105 | } | |
106 | ||
107 | // | |
108 | // Given a stable forwards recurrence relation: | |
109 | // a f_n-1 + b f_n + c f_n+1 = 0 | |
110 | // returns the ratio f_n / f_n+1 | |
111 | // | |
112 | // Note that in most situations where this would be used, we're relying on | |
113 | // pseudo-convergence, as in most cases f_n will not be minimal as N -> -INF | |
114 | // as long as we reach convergence on the continued-fraction before f_n | |
115 | // switches behaviour, we should be fine. | |
116 | // | |
117 | // Recurrence: a functor that returns a tuple of the factors (a,b,c). | |
118 | // factor: Convergence criteria, should be no less than machine epsilon. | |
119 | // max_iter: Maximum iterations to use solving the continued fraction. | |
120 | // | |
121 | template <class Recurrence, class T> | |
122 | T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter) | |
123 | { | |
124 | boost::math::tools::detail::function_ratio_from_backwards_recurrence_fraction<boost::math::tools::detail::recurrence_reverser<Recurrence, T> > f(r); | |
125 | return boost::math::tools::continued_fraction_a(f, factor, max_iter); | |
126 | } | |
127 | ||
128 | ||
129 | ||
130 | // solves usual recurrence relation for homogeneous | |
131 | // difference equation in stable forward direction | |
132 | // a(n)w(n-1) + b(n)w(n) + c(n)w(n+1) = 0 | |
133 | // | |
134 | // Params: | |
135 | // get_coefs: functor returning a tuple, where | |
136 | // get<0>() is a(n); get<1>() is b(n); get<2>() is c(n); | |
137 | // last_index: index N to be found; | |
138 | // first: w(-1); | |
139 | // second: w(0); | |
140 | // | |
141 | template <class NextCoefs, class T> | |
142 | inline T apply_recurrence_relation_forward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0) | |
143 | { | |
144 | BOOST_MATH_STD_USING | |
145 | using boost::math::tuple; | |
146 | using boost::math::get; | |
147 | ||
148 | T third; | |
149 | T a, b, c; | |
150 | ||
151 | for (unsigned k = 0; k < number_of_steps; ++k) | |
152 | { | |
153 | tie(a, b, c) = get_coefs(k); | |
154 | ||
155 | if ((log_scaling) && | |
156 | ((fabs(tools::max_value<T>() * (c / (a * 2048))) < fabs(first)) | |
157 | || (fabs(tools::max_value<T>() * (c / (b * 2048))) < fabs(second)) | |
158 | || (fabs(tools::min_value<T>() * (c * 2048 / a)) > fabs(first)) | |
159 | || (fabs(tools::min_value<T>() * (c * 2048 / b)) > fabs(second)) | |
160 | )) | |
161 | ||
162 | { | |
163 | // Rescale everything: | |
164 | int log_scale = itrunc(log(fabs(second))); | |
165 | T scale = exp(T(-log_scale)); | |
166 | second *= scale; | |
167 | first *= scale; | |
168 | *log_scaling += log_scale; | |
169 | } | |
170 | // scale each part seperately to avoid spurious overflow: | |
171 | third = (a / -c) * first + (b / -c) * second; | |
172 | BOOST_ASSERT((boost::math::isfinite)(third)); | |
173 | ||
174 | ||
175 | swap(first, second); | |
176 | swap(second, third); | |
177 | } | |
178 | ||
179 | if (previous) | |
180 | *previous = first; | |
181 | ||
182 | return second; | |
183 | } | |
184 | ||
185 | // solves usual recurrence relation for homogeneous | |
186 | // difference equation in stable backward direction | |
187 | // a(n)w(n-1) + b(n)w(n) + c(n)w(n+1) = 0 | |
188 | // | |
189 | // Params: | |
190 | // get_coefs: functor returning a tuple, where | |
191 | // get<0>() is a(n); get<1>() is b(n); get<2>() is c(n); | |
192 | // number_of_steps: index N to be found; | |
193 | // first: w(1); | |
194 | // second: w(0); | |
195 | // | |
196 | template <class T, class NextCoefs> | |
197 | inline T apply_recurrence_relation_backward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0) | |
198 | { | |
199 | BOOST_MATH_STD_USING | |
200 | using boost::math::tuple; | |
201 | using boost::math::get; | |
202 | ||
203 | T next; | |
204 | T a, b, c; | |
205 | ||
206 | for (unsigned k = 0; k < number_of_steps; ++k) | |
207 | { | |
208 | tie(a, b, c) = get_coefs(-static_cast<int>(k)); | |
209 | ||
210 | if ((log_scaling) && | |
211 | ( (fabs(tools::max_value<T>() * (a / b) / 2048) < fabs(second)) | |
212 | || (fabs(tools::max_value<T>() * (a / c) / 2048) < fabs(first)) | |
213 | || (fabs(tools::min_value<T>() * (a / b) * 2048) > fabs(second)) | |
214 | || (fabs(tools::min_value<T>() * (a / c) * 2048) > fabs(first)) | |
215 | )) | |
216 | { | |
217 | // Rescale everything: | |
218 | int log_scale = itrunc(log(fabs(second))); | |
219 | T scale = exp(T(-log_scale)); | |
220 | second *= scale; | |
221 | first *= scale; | |
222 | *log_scaling += log_scale; | |
223 | } | |
224 | // scale each part seperately to avoid spurious overflow: | |
225 | next = (b / -a) * second + (c / -a) * first; | |
226 | BOOST_ASSERT((boost::math::isfinite)(next)); | |
227 | ||
228 | swap(first, second); | |
229 | swap(second, next); | |
230 | } | |
231 | ||
232 | if (previous) | |
233 | *previous = first; | |
234 | ||
235 | return second; | |
236 | } | |
237 | ||
238 | template <class Recurrence> | |
239 | struct forward_recurrence_iterator | |
240 | { | |
241 | typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type; | |
242 | ||
243 | forward_recurrence_iterator(const Recurrence& r, value_type f_n_minus_1, value_type f_n) | |
244 | : f_n_minus_1(f_n_minus_1), f_n(f_n), coef(r), k(0) {} | |
245 | ||
246 | forward_recurrence_iterator(const Recurrence& r, value_type f_n) | |
247 | : f_n(f_n), coef(r), k(0) | |
248 | { | |
249 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<boost::math::policies::policy<> >(); | |
250 | f_n_minus_1 = f_n * boost::math::tools::function_ratio_from_forwards_recurrence(detail::recurrence_offsetter<Recurrence>(r, -1), value_type(boost::math::tools::epsilon<value_type>() * 2), max_iter); | |
251 | boost::math::policies::check_series_iterations<value_type>("forward_recurrence_iterator<>::forward_recurrence_iterator", max_iter, boost::math::policies::policy<>()); | |
252 | } | |
253 | ||
254 | forward_recurrence_iterator& operator++() | |
255 | { | |
256 | using std::swap; | |
257 | value_type a, b, c; | |
258 | boost::math::tie(a, b, c) = coef(k); | |
259 | value_type f_n_plus_1 = a * f_n_minus_1 / -c + b * f_n / -c; | |
260 | swap(f_n_minus_1, f_n); | |
261 | swap(f_n, f_n_plus_1); | |
262 | ++k; | |
263 | return *this; | |
264 | } | |
265 | ||
266 | forward_recurrence_iterator operator++(int) | |
267 | { | |
268 | forward_recurrence_iterator t(*this); | |
269 | ++(*this); | |
270 | return t; | |
271 | } | |
272 | ||
273 | value_type operator*() { return f_n; } | |
274 | ||
275 | value_type f_n_minus_1, f_n; | |
276 | Recurrence coef; | |
277 | int k; | |
278 | }; | |
279 | ||
280 | template <class Recurrence> | |
281 | struct backward_recurrence_iterator | |
282 | { | |
283 | typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type; | |
284 | ||
285 | backward_recurrence_iterator(const Recurrence& r, value_type f_n_plus_1, value_type f_n) | |
286 | : f_n_plus_1(f_n_plus_1), f_n(f_n), coef(r), k(0) {} | |
287 | ||
288 | backward_recurrence_iterator(const Recurrence& r, value_type f_n) | |
289 | : f_n(f_n), coef(r), k(0) | |
290 | { | |
291 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<boost::math::policies::policy<> >(); | |
292 | f_n_plus_1 = f_n * boost::math::tools::function_ratio_from_backwards_recurrence(detail::recurrence_offsetter<Recurrence>(r, 1), value_type(boost::math::tools::epsilon<value_type>() * 2), max_iter); | |
293 | boost::math::policies::check_series_iterations<value_type>("backward_recurrence_iterator<>::backward_recurrence_iterator", max_iter, boost::math::policies::policy<>()); | |
294 | } | |
295 | ||
296 | backward_recurrence_iterator& operator++() | |
297 | { | |
298 | using std::swap; | |
299 | value_type a, b, c; | |
300 | boost::math::tie(a, b, c) = coef(k); | |
301 | value_type f_n_minus_1 = c * f_n_plus_1 / -a + b * f_n / -a; | |
302 | swap(f_n_plus_1, f_n); | |
303 | swap(f_n, f_n_minus_1); | |
304 | --k; | |
305 | return *this; | |
306 | } | |
307 | ||
308 | backward_recurrence_iterator operator++(int) | |
309 | { | |
310 | backward_recurrence_iterator t(*this); | |
311 | ++(*this); | |
312 | return t; | |
313 | } | |
314 | ||
315 | value_type operator*() { return f_n; } | |
316 | ||
317 | value_type f_n_plus_1, f_n; | |
318 | Recurrence coef; | |
319 | int k; | |
320 | }; | |
321 | ||
322 | } | |
323 | } | |
324 | } // namespaces | |
325 | ||
326 | #endif // BOOST_MATH_TOOLS_RECURRENCE_HPP_ |