boost/math/tools/recurrence.hpp

   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
  14 #ifdef BOOST_NO_CXX11_HDR_TUPLE
  15 #error "This header requires C++11 support"
  16 #endif
  17
  18
  19 namespace boost {
  20    namespace math {
  21       namespace tools {
  22          namespace detail{
  23
  24             //
  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.
  28             // and:
  29             // COMPUTATIONAL ASPECTS OF THREE-TERM RECURRENCE RELATIONS
  30             // WALTER GAUTSCHI
  31             // SIAM REVIEW Vol. 9, No. 1, January, 1967
  32             //
  33             template <class Recurrence>
  34             struct function_ratio_from_backwards_recurrence_fraction
  35             {
  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) {}
  39
  40                result_type operator()()
  41                {
  42                   value_type a, b, c;
  43                   boost::math::tie(a, b, c) = r(k);
  44                   ++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);
  50                }
  51
  52             private:
  53                function_ratio_from_backwards_recurrence_fraction operator=(const function_ratio_from_backwards_recurrence_fraction&);
  54
  55                Recurrence r;
  56                int k;
  57             };
  58
  59             template <class R, class T>
  60             struct recurrence_reverser
  61             {
  62                recurrence_reverser(const R& r) : r(r) {}
  63                boost::math::tuple<T, T, T> operator()(int i)
  64                {
  65                   using std::swap;
  66                   boost::math::tuple<T, T, T> t = r(-i);
  67                   swap(boost::math::get<0>(t), boost::math::get<2>(t));
  68                   return t;
  69                }
  70                R r;
  71             };
  72
  73             template <class Recurrence>
  74             struct recurrence_offsetter
  75             {
  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)
  79                {
  80                   return r(i + k);
  81                }
  82             private:
  83                Recurrence r;
  84                int k;
  85             };
  86
  87
  88
  89          }  // namespace detail
  90
  91          //
  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
  95          //
  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.
  99          //
 100          template <class Recurrence, class T>
 101          T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter)
 102          {
 103             detail::function_ratio_from_backwards_recurrence_fraction<Recurrence> f(r);
 104             return boost::math::tools::continued_fraction_a(f, factor, max_iter);
 105          }
 106
 107          //
 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
 111          //
 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.
 116          //
 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.
 120          //
 121          template <class Recurrence, class T>
 122          T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter)
 123          {
 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);
 126          }
 127
 128
 129
 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
 133          //
 134          // Params:
 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;
 138          // first: w(-1);
 139          // second: w(0);
 140          //
 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)
 143          {
 144             BOOST_MATH_STD_USING
 145             using boost::math::tuple;
 146             using boost::math::get;
 147
 148             T third;
 149             T a, b, c;
 150
 151             for (unsigned k = 0; k < number_of_steps; ++k)
 152             {
 153                tie(a, b, c) = get_coefs(k);
 154
 155                if ((log_scaling) &&
 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))
 160                      ))
 161
 162                {
 163                   // Rescale everything:
 164                   int log_scale = itrunc(log(fabs(second)));
 165                   T scale = exp(T(-log_scale));
 166                   second *= scale;
 167                   first *= scale;
 168                   *log_scaling += log_scale;
 169                }
 170                // scale each part seperately to avoid spurious overflow:
 171                third = (a / -c) * first + (b / -c) * second;
 172                BOOST_ASSERT((boost::math::isfinite)(third));
 173
 174
 175                swap(first, second);
 176                swap(second, third);
 177             }
 178
 179             if (previous)
 180                *previous = first;
 181
 182             return second;
 183          }
 184
 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
 188          //
 189          // Params:
 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;
 193          // first: w(1);
 194          // second: w(0);
 195          //
 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)
 198          {
 199             BOOST_MATH_STD_USING
 200             using boost::math::tuple;
 201             using boost::math::get;
 202
 203             T next;
 204             T a, b, c;
 205
 206             for (unsigned k = 0; k < number_of_steps; ++k)
 207             {
 208                tie(a, b, c) = get_coefs(-static_cast<int>(k));
 209
 210                if ((log_scaling) &&
 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))
 215                   ))
 216                {
 217                   // Rescale everything:
 218                   int log_scale = itrunc(log(fabs(second)));
 219                   T scale = exp(T(-log_scale));
 220                   second *= scale;
 221                   first *= scale;
 222                   *log_scaling += log_scale;
 223                }
 224                // scale each part seperately to avoid spurious overflow:
 225                next = (b / -a) * second + (c / -a) * first;
 226                BOOST_ASSERT((boost::math::isfinite)(next));
 227
 228                swap(first, second);
 229                swap(second, next);
 230             }
 231
 232             if (previous)
 233                *previous = first;
 234
 235             return second;
 236          }
 237
 238          template <class Recurrence>
 239          struct forward_recurrence_iterator
 240          {
 241             typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
 242
 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) {}
 245
 246             forward_recurrence_iterator(const Recurrence& r, value_type f_n)
 247                : f_n(f_n), coef(r), k(0)
 248             {
 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<>());
 252             }
 253
 254             forward_recurrence_iterator& operator++()
 255             {
 256                using std::swap;
 257                value_type a, b, c;
 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);
 262                ++k;
 263                return *this;
 264             }
 265
 266             forward_recurrence_iterator operator++(int)
 267             {
 268                forward_recurrence_iterator t(*this);
 269                ++(*this);
 270                return t;
 271             }
 272
 273             value_type operator*() { return f_n; }
 274
 275             value_type f_n_minus_1, f_n;
 276             Recurrence coef;
 277             int k;
 278          };
 279
 280          template <class Recurrence>
 281          struct backward_recurrence_iterator
 282          {
 283             typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
 284
 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) {}
 287
 288             backward_recurrence_iterator(const Recurrence& r, value_type f_n)
 289                : f_n(f_n), coef(r), k(0)
 290             {
 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<>());
 294             }
 295
 296             backward_recurrence_iterator& operator++()
 297             {
 298                using std::swap;
 299                value_type a, b, c;
 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);
 304                --k;
 305                return *this;
 306             }
 307
 308             backward_recurrence_iterator operator++(int)
 309             {
 310                backward_recurrence_iterator t(*this);
 311                ++(*this);
 312                return t;
 313             }
 314
 315             value_type operator*() { return f_n; }
 316
 317             value_type f_n_plus_1, f_n;
 318             Recurrence coef;
 319             int k;
 320          };
 321
 322       }
 323    }
 324 } // namespaces
 325
 326 #endif // BOOST_MATH_TOOLS_RECURRENCE_HPP_