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)
6 #ifndef BOOST_MATH_TOOLS_RECURRENCE_HPP_
7 #define BOOST_MATH_TOOLS_RECURRENCE_HPP_
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>
14 #ifdef BOOST_NO_CXX11_HDR_TUPLE
15 #error "This header requires C++11 support"
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.
29 // COMPUTATIONAL ASPECTS OF THREE-TERM RECURRENCE RELATIONS
31 // SIAM REVIEW Vol. 9, No. 1, January, 1967
33 template <class Recurrence>
34 struct function_ratio_from_backwards_recurrence_fraction
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) {}
40 result_type operator()()
43 boost::math::tie(a, b, c) = r(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);
53 function_ratio_from_backwards_recurrence_fraction operator=(const function_ratio_from_backwards_recurrence_fraction&);
59 template <class R, class T>
60 struct recurrence_reverser
62 recurrence_reverser(const R& r) : r(r) {}
63 boost::math::tuple<T, T, T> operator()(int i)
66 boost::math::tuple<T, T, T> t = r(-i);
67 swap(boost::math::get<0>(t), boost::math::get<2>(t));
73 template <class Recurrence>
74 struct recurrence_offsetter
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)
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
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.
100 template <class Recurrence, class T>
101 T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter)
103 detail::function_ratio_from_backwards_recurrence_fraction<Recurrence> f(r);
104 return boost::math::tools::continued_fraction_a(f, factor, max_iter);
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
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.
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.
121 template <class Recurrence, class T>
122 T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter)
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);
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
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;
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)
145 using boost::math::tuple;
146 using boost::math::get;
151 for (unsigned k = 0; k < number_of_steps; ++k)
153 tie(a, b, c) = get_coefs(k);
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))
163 // Rescale everything:
164 int log_scale = itrunc(log(fabs(second)));
165 T scale = exp(T(-log_scale));
168 *log_scaling += log_scale;
170 // scale each part seperately to avoid spurious overflow:
171 third = (a / -c) * first + (b / -c) * second;
172 BOOST_ASSERT((boost::math::isfinite)(third));
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
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;
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)
200 using boost::math::tuple;
201 using boost::math::get;
206 for (unsigned k = 0; k < number_of_steps; ++k)
208 tie(a, b, c) = get_coefs(-static_cast<int>(k));
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))
217 // Rescale everything:
218 int log_scale = itrunc(log(fabs(second)));
219 T scale = exp(T(-log_scale));
222 *log_scaling += log_scale;
224 // scale each part seperately to avoid spurious overflow:
225 next = (b / -a) * second + (c / -a) * first;
226 BOOST_ASSERT((boost::math::isfinite)(next));
238 template <class Recurrence>
239 struct forward_recurrence_iterator
241 typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
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) {}
246 forward_recurrence_iterator(const Recurrence& r, value_type f_n)
247 : f_n(f_n), coef(r), k(0)
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<>());
254 forward_recurrence_iterator& operator++()
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);
266 forward_recurrence_iterator operator++(int)
268 forward_recurrence_iterator t(*this);
273 value_type operator*() { return f_n; }
275 value_type f_n_minus_1, f_n;
280 template <class Recurrence>
281 struct backward_recurrence_iterator
283 typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
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) {}
288 backward_recurrence_iterator(const Recurrence& r, value_type f_n)
289 : f_n(f_n), coef(r), k(0)
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<>());
296 backward_recurrence_iterator& operator++()
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);
308 backward_recurrence_iterator operator++(int)
310 backward_recurrence_iterator t(*this);
315 value_type operator*() { return f_n; }
317 value_type f_n_plus_1, f_n;
326 #endif // BOOST_MATH_TOOLS_RECURRENCE_HPP_