math/tools/series.hpp

0001 //  (C) Copyright John Maddock 2005-2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED
0007 #define BOOST_MATH_TOOLS_SERIES_INCLUDED
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <cmath>
0014 #include <cstdint>
0015 #include <limits>
0016 #include <boost/math/tools/config.hpp>
0017
0018 namespace boost{ namespace math{ namespace tools{
0019
0020 //
0021 // Simple series summation come first:
0022 //
0023 template <class Functor, class U, class V>
0024 inline typename Functor::result_type sum_series(Functor& func, const U& factor, std::uintmax_t& max_terms, const V& init_value) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0025 {
0026    BOOST_MATH_STD_USING
0027
0028    typedef typename Functor::result_type result_type;
0029
0030    std::uintmax_t counter = max_terms;
0031
0032    result_type result = init_value;
0033    result_type next_term;
0034    do{
0035       next_term = func();
0036       result += next_term;
0037    }
0038    while((abs(factor * result) < abs(next_term)) && --counter);
0039
0040    // set max_terms to the actual number of terms of the series evaluated:
0041    max_terms = max_terms - counter;
0042
0043    return result;
0044 }
0045
0046 template <class Functor, class U>
0047 inline typename Functor::result_type sum_series(Functor& func, const U& factor, std::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0048 {
0049    typename Functor::result_type init_value = 0;
0050    return sum_series(func, factor, max_terms, init_value);
0051 }
0052
0053 template <class Functor, class U>
0054 inline typename Functor::result_type sum_series(Functor& func, int bits, std::uintmax_t& max_terms, const U& init_value) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0055 {
0056    BOOST_MATH_STD_USING
0057    typedef typename Functor::result_type result_type;
0058    result_type factor = ldexp(result_type(1), 1 - bits);
0059    return sum_series(func, factor, max_terms, init_value);
0060 }
0061
0062 template <class Functor>
0063 inline typename Functor::result_type sum_series(Functor& func, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0064 {
0065    BOOST_MATH_STD_USING
0066    typedef typename Functor::result_type result_type;
0067    std::uintmax_t iters = (std::numeric_limits<std::uintmax_t>::max)();
0068    result_type init_val = 0;
0069    return sum_series(func, bits, iters, init_val);
0070 }
0071
0072 template <class Functor>
0073 inline typename Functor::result_type sum_series(Functor& func, int bits, std::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0074 {
0075    BOOST_MATH_STD_USING
0076    typedef typename Functor::result_type result_type;
0077    result_type init_val = 0;
0078    return sum_series(func, bits, max_terms, init_val);
0079 }
0080
0081 template <class Functor, class U>
0082 inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0083 {
0084    BOOST_MATH_STD_USING
0085    std::uintmax_t iters = (std::numeric_limits<std::uintmax_t>::max)();
0086    return sum_series(func, bits, iters, init_value);
0087 }
0088 //
0089 // Checked summation:
0090 //
0091 template <class Functor, class U, class V>
0092 inline typename Functor::result_type checked_sum_series(Functor& func, const U& factor, std::uintmax_t& max_terms, const V& init_value, V& norm) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0093 {
0094    BOOST_MATH_STD_USING
0095
0096    typedef typename Functor::result_type result_type;
0097
0098    std::uintmax_t counter = max_terms;
0099
0100    result_type result = init_value;
0101    result_type next_term;
0102    do {
0103       next_term = func();
0104       result += next_term;
0105       norm += fabs(next_term);
0106    } while ((abs(factor * result) < abs(next_term)) && --counter);
0107
0108    // set max_terms to the actual number of terms of the series evaluated:
0109    max_terms = max_terms - counter;
0110
0111    return result;
0112 }
0113
0114
0115 //
0116 // Algorithm kahan_sum_series invokes Functor func until the N'th
0117 // term is too small to have any effect on the total, the terms
0118 // are added using the Kahan summation method.
0119 //
0120 // CAUTION: Optimizing compilers combined with extended-precision
0121 // machine registers conspire to render this algorithm partly broken:
0122 // double rounding of intermediate terms (first to a long double machine
0123 // register, and then to a double result) cause the rounding error computed
0124 // by the algorithm to be off by up to 1ulp.  However this occurs rarely, and
0125 // in any case the result is still much better than a naive summation.
0126 //
0127 template <class Functor>
0128 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0129 {
0130    BOOST_MATH_STD_USING
0131
0132    typedef typename Functor::result_type result_type;
0133
0134    result_type factor = pow(result_type(2), result_type(bits));
0135    result_type result = func();
0136    result_type next_term, y, t;
0137    result_type carry = 0;
0138    do{
0139       next_term = func();
0140       y = next_term - carry;
0141       t = result + y;
0142       carry = t - result;
0143       carry -= y;
0144       result = t;
0145    }
0146    while(fabs(result) < fabs(factor * next_term));
0147    return result;
0148 }
0149
0150 template <class Functor>
0151 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, std::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
0152 {
0153    BOOST_MATH_STD_USING
0154
0155    typedef typename Functor::result_type result_type;
0156
0157    std::uintmax_t counter = max_terms;
0158
0159    result_type factor = ldexp(result_type(1), bits);
0160    result_type result = func();
0161    result_type next_term, y, t;
0162    result_type carry = 0;
0163    do{
0164       next_term = func();
0165       y = next_term - carry;
0166       t = result + y;
0167       carry = t - result;
0168       carry -= y;
0169       result = t;
0170    }
0171    while((fabs(result) < fabs(factor * next_term)) && --counter);
0172
0173    // set max_terms to the actual number of terms of the series evaluated:
0174    max_terms = max_terms - counter;
0175
0176    return result;
0177 }
0178
0179 } // namespace tools
0180 } // namespace math
0181 } // namespace boost
0182
0183 #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED
0184