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