boost/math/tools/series.hpp

   1 //  (C) Copyright John Maddock 2005-2006.
   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_SERIES_INCLUDED
   7 #define BOOST_MATH_TOOLS_SERIES_INCLUDED
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/config/no_tr1/cmath.hpp>
  14 #include <boost/cstdint.hpp>
  15 #include <boost/limits.hpp>
  16 #include <boost/math/tools/config.hpp>
  17
  18 namespace boost{ namespace math{ namespace tools{
  19
  20 //
  21 // Simple series summation come first:
  22 //
  23 template <class Functor, class U, class V>
  24 inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  25 {
  26    BOOST_MATH_STD_USING
  27
  28    typedef typename Functor::result_type result_type;
  29
  30    boost::uintmax_t counter = max_terms;
  31
  32    result_type result = init_value;
  33    result_type next_term;
  34    do{
  35       next_term = func();
  36       result += next_term;
  37    }
  38    while((abs(factor * result) < abs(next_term)) && --counter);
  39
  40    // set max_terms to the actual number of terms of the series evaluated:
  41    max_terms = max_terms - counter;
  42
  43    return result;
  44 }
  45
  46 template <class Functor, class U>
  47 inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  48 {
  49    typename Functor::result_type init_value = 0;
  50    return sum_series(func, factor, max_terms, init_value);
  51 }
  52
  53 template <class Functor, class U>
  54 inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  55 {
  56    BOOST_MATH_STD_USING
  57    typedef typename Functor::result_type result_type;
  58    result_type factor = ldexp(result_type(1), 1 - bits);
  59    return sum_series(func, factor, max_terms, init_value);
  60 }
  61
  62 template <class Functor>
  63 inline typename Functor::result_type sum_series(Functor& func, int bits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  64 {
  65    BOOST_MATH_STD_USING
  66    typedef typename Functor::result_type result_type;
  67    boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
  68    result_type init_val = 0;
  69    return sum_series(func, bits, iters, init_val);
  70 }
  71
  72 template <class Functor>
  73 inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  74 {
  75    BOOST_MATH_STD_USING
  76    typedef typename Functor::result_type result_type;
  77    result_type init_val = 0;
  78    return sum_series(func, bits, max_terms, init_val);
  79 }
  80
  81 template <class Functor, class U>
  82 inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  83 {
  84    BOOST_MATH_STD_USING
  85    boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
  86    return sum_series(func, bits, iters, init_value);
  87 }
  88 //
  89 // Checked summation:
  90 //
  91 template <class Functor, class U, class V>
  92 inline typename Functor::result_type checked_sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value, V& norm) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  93 {
  94    BOOST_MATH_STD_USING
  95
  96    typedef typename Functor::result_type result_type;
  97
  98    boost::uintmax_t counter = max_terms;
  99
 100    result_type result = init_value;
 101    result_type next_term;
 102    do {
 103       next_term = func();
 104       result += next_term;
 105       norm += fabs(next_term);
 106    } while ((abs(factor * result) < abs(next_term)) && --counter);
 107
 108    // set max_terms to the actual number of terms of the series evaluated:
 109    max_terms = max_terms - counter;
 110
 111    return result;
 112 }
 113
 114
 115 //
 116 // Algorithm kahan_sum_series invokes Functor func until the N'th
 117 // term is too small to have any effect on the total, the terms
 118 // are added using the Kahan summation method.
 119 //
 120 // CAUTION: Optimizing compilers combined with extended-precision
 121 // machine registers conspire to render this algorithm partly broken:
 122 // double rounding of intermediate terms (first to a long double machine
 123 // register, and then to a double result) cause the rounding error computed
 124 // by the algorithm to be off by up to 1ulp.  However this occurs rarely, and
 125 // in any case the result is still much better than a naive summation.
 126 //
 127 template <class Functor>
 128 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
 129 {
 130    BOOST_MATH_STD_USING
 131
 132    typedef typename Functor::result_type result_type;
 133
 134    result_type factor = pow(result_type(2), bits);
 135    result_type result = func();
 136    result_type next_term, y, t;
 137    result_type carry = 0;
 138    do{
 139       next_term = func();
 140       y = next_term - carry;
 141       t = result + y;
 142       carry = t - result;
 143       carry -= y;
 144       result = t;
 145    }
 146    while(fabs(result) < fabs(factor * next_term));
 147    return result;
 148 }
 149
 150 template <class Functor>
 151 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
 152 {
 153    BOOST_MATH_STD_USING
 154
 155    typedef typename Functor::result_type result_type;
 156
 157    boost::uintmax_t counter = max_terms;
 158
 159    result_type factor = ldexp(result_type(1), bits);
 160    result_type result = func();
 161    result_type next_term, y, t;
 162    result_type carry = 0;
 163    do{
 164       next_term = func();
 165       y = next_term - carry;
 166       t = result + y;
 167       carry = t - result;
 168       carry -= y;
 169       result = t;
 170    }
 171    while((fabs(result) < fabs(factor * next_term)) && --counter);
 172
 173    // set max_terms to the actual number of terms of the series evaluated:
 174    max_terms = max_terms - counter;
 175
 176    return result;
 177 }
 178
 179 } // namespace tools
 180 } // namespace math
 181 } // namespace boost
 182
 183 #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED
 184