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