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)
  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((fabs(factor * result) < fabs(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)
  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)
  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)
  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)
  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)
  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 //
  90 // Algorithm kahan_sum_series invokes Functor func until the N'th
  91 // term is too small to have any effect on the total, the terms
  92 // are added using the Kahan summation method.
  93 //
  94 // CAUTION: Optimizing compilers combined with extended-precision
  95 // machine registers conspire to render this algorithm partly broken:
  96 // double rounding of intermediate terms (first to a long double machine
  97 // register, and then to a double result) cause the rounding error computed
  98 // by the algorithm to be off by up to 1ulp.  However this occurs rarely, and
  99 // in any case the result is still much better than a naive summation.
 100 //
 101 template <class Functor>
 102 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits)
 103 {
 104    BOOST_MATH_STD_USING
 105
 106    typedef typename Functor::result_type result_type;
 107
 108    result_type factor = pow(result_type(2), bits);
 109    result_type result = func();
 110    result_type next_term, y, t;
 111    result_type carry = 0;
 112    do{
 113       next_term = func();
 114       y = next_term - carry;
 115       t = result + y;
 116       carry = t - result;
 117       carry -= y;
 118       result = t;
 119    }
 120    while(fabs(result) < fabs(factor * next_term));
 121    return result;
 122 }
 123
 124 template <class Functor>
 125 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
 126 {
 127    BOOST_MATH_STD_USING
 128
 129    typedef typename Functor::result_type result_type;
 130
 131    boost::uintmax_t counter = max_terms;
 132
 133    result_type factor = ldexp(result_type(1), bits);
 134    result_type result = func();
 135    result_type next_term, y, t;
 136    result_type carry = 0;
 137    do{
 138       next_term = func();
 139       y = next_term - carry;
 140       t = result + y;
 141       carry = t - result;
 142       carry -= y;
 143       result = t;
 144    }
 145    while((fabs(result) < fabs(factor * next_term)) && --counter);
 146
 147    // set max_terms to the actual number of terms of the series evaluated:
 148    max_terms = max_terms - counter;
 149
 150    return result;
 151 }
 152
 153 } // namespace tools
 154 } // namespace math
 155 } // namespace boost
 156
 157 #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED
 158