math/tools/recurrence.hpp

0001 //  (C) Copyright Anton Bikineev 2014
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_RECURRENCE_HPP_
0007 #define BOOST_MATH_TOOLS_RECURRENCE_HPP_
0008
0009 #include <type_traits>
0010 #include <tuple>
0011 #include <utility>
0012 #include <boost/math/tools/config.hpp>
0013 #include <boost/math/tools/precision.hpp>
0014 #include <boost/math/tools/tuple.hpp>
0015 #include <boost/math/tools/fraction.hpp>
0016 #include <boost/math/tools/cxx03_warn.hpp>
0017 #include <boost/math/tools/assert.hpp>
0018
0019 namespace boost {
0020    namespace math {
0021       namespace tools {
0022          namespace detail{
0023
0024             //
0025             // Function ratios directly from recurrence relations:
0026             // H. Shintan, Note on Miller's recurrence algorithm, J. Sci. Hiroshima Univ. Ser. A-I
0027             // Math., 29 (1965), pp. 121 - 133.
0028             // and:
0029             // COMPUTATIONAL ASPECTS OF THREE-TERM RECURRENCE RELATIONS
0030             // WALTER GAUTSCHI
0031             // SIAM REVIEW Vol. 9, No. 1, January, 1967
0032             //
0033             template <class Recurrence>
0034             struct function_ratio_from_backwards_recurrence_fraction
0035             {
0036                typedef typename std::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
0037                typedef std::pair<value_type, value_type> result_type;
0038                function_ratio_from_backwards_recurrence_fraction(const Recurrence& r) : r(r), k(0) {}
0039
0040                result_type operator()()
0041                {
0042                   value_type a, b, c;
0043                   std::tie(a, b, c) = r(k);
0044                   ++k;
0045                   // an and bn defined as per Gauchi 1.16, not the same
0046                   // as the usual continued fraction a' and b's.
0047                   value_type bn = a / c;
0048                   value_type an = b / c;
0049                   return result_type(-bn, an);
0050                }
0051
0052             private:
0053                function_ratio_from_backwards_recurrence_fraction operator=(const function_ratio_from_backwards_recurrence_fraction&) = delete;
0054
0055                Recurrence r;
0056                int k;
0057             };
0058
0059             template <class R, class T>
0060             struct recurrence_reverser
0061             {
0062                recurrence_reverser(const R& r) : r(r) {}
0063                std::tuple<T, T, T> operator()(int i)
0064                {
0065                   using std::swap;
0066                   std::tuple<T, T, T> t = r(-i);
0067                   swap(std::get<0>(t), std::get<2>(t));
0068                   return t;
0069                }
0070                R r;
0071             };
0072
0073             template <class Recurrence>
0074             struct recurrence_offsetter
0075             {
0076                typedef decltype(std::declval<Recurrence&>()(0)) result_type;
0077                recurrence_offsetter(Recurrence const& rr, int offset) : r(rr), k(offset) {}
0078                result_type operator()(int i)
0079                {
0080                   return r(i + k);
0081                }
0082             private:
0083                Recurrence r;
0084                int k;
0085             };
0086
0087
0088
0089          }  // namespace detail
0090
0091          //
0092          // Given a stable backwards recurrence relation:
0093          // a f_n-1 + b f_n + c f_n+1 = 0
0094          // returns the ratio f_n / f_n-1
0095          //
0096          // Recurrence: a functor that returns a tuple of the factors (a,b,c).
0097          // factor:     Convergence criteria, should be no less than machine epsilon.
0098          // max_iter:   Maximum iterations to use solving the continued fraction.
0099          //
0100          template <class Recurrence, class T>
0101          T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, std::uintmax_t& max_iter)
0102          {
0103             detail::function_ratio_from_backwards_recurrence_fraction<Recurrence> f(r);
0104             return boost::math::tools::continued_fraction_a(f, factor, max_iter);
0105          }
0106
0107          //
0108          // Given a stable forwards recurrence relation:
0109          // a f_n-1 + b f_n + c f_n+1 = 0
0110          // returns the ratio f_n / f_n+1
0111          //
0112          // Note that in most situations where this would be used, we're relying on
0113          // pseudo-convergence, as in most cases f_n will not be minimal as N -> -INF
0114          // as long as we reach convergence on the continued-fraction before f_n
0115          // switches behaviour, we should be fine.
0116          //
0117          // Recurrence: a functor that returns a tuple of the factors (a,b,c).
0118          // factor:     Convergence criteria, should be no less than machine epsilon.
0119          // max_iter:   Maximum iterations to use solving the continued fraction.
0120          //
0121          template <class Recurrence, class T>
0122          T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, std::uintmax_t& max_iter)
0123          {
0124             boost::math::tools::detail::function_ratio_from_backwards_recurrence_fraction<boost::math::tools::detail::recurrence_reverser<Recurrence, T> > f(r);
0125             return boost::math::tools::continued_fraction_a(f, factor, max_iter);
0126          }
0127
0128
0129
0130          // solves usual recurrence relation for homogeneous
0131          // difference equation in stable forward direction
0132          // a(n)w(n-1) + b(n)w(n) + c(n)w(n+1) = 0
0133          //
0134          // Params:
0135          // get_coefs: functor returning a tuple, where
0136          //            get<0>() is a(n); get<1>() is b(n); get<2>() is c(n);
0137          // last_index: index N to be found;
0138          // first: w(-1);
0139          // second: w(0);
0140          //
0141          template <class NextCoefs, class T>
0142          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)
0143          {
0144             BOOST_MATH_STD_USING
0145             using std::tuple;
0146             using std::get;
0147             using std::swap;
0148
0149             T third;
0150             T a, b, c;
0151
0152             for (unsigned k = 0; k < number_of_steps; ++k)
0153             {
0154                tie(a, b, c) = get_coefs(k);
0155
0156                if ((log_scaling) &&
0157                   ((fabs(tools::max_value<T>() * (c / (a * 2048))) < fabs(first))
0158                      || (fabs(tools::max_value<T>() * (c / (b * 2048))) < fabs(second))
0159                      || (fabs(tools::min_value<T>() * (c * 2048 / a)) > fabs(first))
0160                      || (fabs(tools::min_value<T>() * (c * 2048 / b)) > fabs(second))
0161                      ))
0162
0163                {
0164                   // Rescale everything:
0165                   long long log_scale = lltrunc(log(fabs(second)));
0166                   T scale = exp(T(-log_scale));
0167                   second *= scale;
0168                   first *= scale;
0169                   *log_scaling += log_scale;
0170                }
0171                // scale each part separately to avoid spurious overflow:
0172                third = (a / -c) * first + (b / -c) * second;
0173                BOOST_MATH_ASSERT((boost::math::isfinite)(third));
0174
0175
0176                swap(first, second);
0177                swap(second, third);
0178             }
0179
0180             if (previous)
0181                *previous = first;
0182
0183             return second;
0184          }
0185
0186          // solves usual recurrence relation for homogeneous
0187          // difference equation in stable backward direction
0188          // a(n)w(n-1) + b(n)w(n) + c(n)w(n+1) = 0
0189          //
0190          // Params:
0191          // get_coefs: functor returning a tuple, where
0192          //            get<0>() is a(n); get<1>() is b(n); get<2>() is c(n);
0193          // number_of_steps: index N to be found;
0194          // first: w(1);
0195          // second: w(0);
0196          //
0197          template <class T, class NextCoefs>
0198          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)
0199          {
0200             BOOST_MATH_STD_USING
0201             using std::tuple;
0202             using std::get;
0203             using std::swap;
0204
0205             T next;
0206             T a, b, c;
0207
0208             for (unsigned k = 0; k < number_of_steps; ++k)
0209             {
0210                tie(a, b, c) = get_coefs(-static_cast<int>(k));
0211
0212                if ((log_scaling) && (second != 0) &&
0213                   ( (fabs(tools::max_value<T>() * (a / b) / 2048) < fabs(second))
0214                      || (fabs(tools::max_value<T>() * (a / c) / 2048) < fabs(first))
0215                      || (fabs(tools::min_value<T>() * (a / b) * 2048) > fabs(second))
0216                      || (fabs(tools::min_value<T>() * (a / c) * 2048) > fabs(first))
0217                   ))
0218                {
0219                   // Rescale everything:
0220                   int log_scale = itrunc(log(fabs(second)));
0221                   T scale = exp(T(-log_scale));
0222                   second *= scale;
0223                   first *= scale;
0224                   *log_scaling += log_scale;
0225                }
0226                // scale each part separately to avoid spurious overflow:
0227                next = (b / -a) * second + (c / -a) * first;
0228                BOOST_MATH_ASSERT((boost::math::isfinite)(next));
0229
0230                swap(first, second);
0231                swap(second, next);
0232             }
0233
0234             if (previous)
0235                *previous = first;
0236
0237             return second;
0238          }
0239
0240          template <class Recurrence>
0241          struct forward_recurrence_iterator
0242          {
0243             typedef typename std::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
0244
0245             forward_recurrence_iterator(const Recurrence& r, value_type f_n_minus_1, value_type f_n)
0246                : f_n_minus_1(f_n_minus_1), f_n(f_n), coef(r), k(0) {}
0247
0248             forward_recurrence_iterator(const Recurrence& r, value_type f_n)
0249                : f_n(f_n), coef(r), k(0)
0250             {
0251                std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<boost::math::policies::policy<> >();
0252                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);
0253                boost::math::policies::check_series_iterations<value_type>("forward_recurrence_iterator<>::forward_recurrence_iterator", max_iter, boost::math::policies::policy<>());
0254             }
0255
0256             forward_recurrence_iterator& operator++()
0257             {
0258                using std::swap;
0259                value_type a, b, c;
0260                std::tie(a, b, c) = coef(k);
0261                value_type f_n_plus_1 = a * f_n_minus_1 / -c + b * f_n / -c;
0262                swap(f_n_minus_1, f_n);
0263                swap(f_n, f_n_plus_1);
0264                ++k;
0265                return *this;
0266             }
0267
0268             forward_recurrence_iterator operator++(int)
0269             {
0270                forward_recurrence_iterator t(*this);
0271                ++(*this);
0272                return t;
0273             }
0274
0275             value_type operator*() { return f_n; }
0276
0277             value_type f_n_minus_1, f_n;
0278             Recurrence coef;
0279             int k;
0280          };
0281
0282          template <class Recurrence>
0283          struct backward_recurrence_iterator
0284          {
0285             typedef typename std::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
0286
0287             backward_recurrence_iterator(const Recurrence& r, value_type f_n_plus_1, value_type f_n)
0288                : f_n_plus_1(f_n_plus_1), f_n(f_n), coef(r), k(0) {}
0289
0290             backward_recurrence_iterator(const Recurrence& r, value_type f_n)
0291                : f_n(f_n), coef(r), k(0)
0292             {
0293                std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<boost::math::policies::policy<> >();
0294                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);
0295                boost::math::policies::check_series_iterations<value_type>("backward_recurrence_iterator<>::backward_recurrence_iterator", max_iter, boost::math::policies::policy<>());
0296             }
0297
0298             backward_recurrence_iterator& operator++()
0299             {
0300                using std::swap;
0301                value_type a, b, c;
0302                std::tie(a, b, c) = coef(k);
0303                value_type f_n_minus_1 = c * f_n_plus_1 / -a + b * f_n / -a;
0304                swap(f_n_plus_1, f_n);
0305                swap(f_n, f_n_minus_1);
0306                --k;
0307                return *this;
0308             }
0309
0310             backward_recurrence_iterator operator++(int)
0311             {
0312                backward_recurrence_iterator t(*this);
0313                ++(*this);
0314                return t;
0315             }
0316
0317             value_type operator*() { return f_n; }
0318
0319             value_type f_n_plus_1, f_n;
0320             Recurrence coef;
0321             int k;
0322          };
0323
0324       }
0325    }
0326 } // namespaces
0327
0328 #endif // BOOST_MATH_TOOLS_RECURRENCE_HPP_