ceph/src/boost/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 #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;
  39                function_ratio_from_backwards_recurrence_fraction(const Recurrence& r) : r(r), k(0) {}
  40
  41                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             {
  63                recurrence_reverser(const R& r) : r(r) {}
  64                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;
  78                recurrence_offsetter(Recurrence const& rr, int offset) : r(rr), k(offset) {}
  79                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>
 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>
 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>
 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>
 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
 244             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
 247             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
 255             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
 267             forward_recurrence_iterator operator++(int)
 268             {
 269                forward_recurrence_iterator t(*this);
 270                ++(*this);
 271                return t;
 272             }
 273
 274             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
 286             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
 289             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
 297             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
 309             backward_recurrence_iterator operator++(int)
 310             {
 311                backward_recurrence_iterator t(*this);
 312                ++(*this);
 313                return t;
 314             }
 315
 316             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_