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