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