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
0014 #include <boost/math/tools/config.hpp>
0015 #include <boost/math/tools/numeric_limits.hpp>
0016 #include <boost/math/tools/cstdint.hpp>
0017 #include <boost/math/tools/type_traits.hpp>
0018
0019 namespace boost{ namespace math{ namespace tools{
0020
0021 //
0022 // Simple series summation come first:
0023 //
0024 template <class Functor, class U, class V>
0025 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::math::uintmax_t& max_terms, const V& init_value) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0026 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0027 && noexcept(std::declval<Functor>()())
0028 #endif
0029 )
0030 {
0031    BOOST_MATH_STD_USING
0032
0033    typedef typename Functor::result_type result_type;
0034
0035    boost::math::uintmax_t counter = max_terms;
0036
0037    result_type result = init_value;
0038    result_type next_term;
0039    do{
0040       next_term = func();
0041       result += next_term;
0042    }
0043    while((abs(factor * result) < abs(next_term)) && --counter);
0044
0045    // set max_terms to the actual number of terms of the series evaluated:
0046    max_terms = max_terms - counter;
0047
0048    return result;
0049 }
0050
0051 template <class Functor, class U>
0052 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0053 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0054 && noexcept(std::declval<Functor>()())
0055 #endif
0056 )
0057 {
0058    typename Functor::result_type init_value = 0;
0059    return sum_series(func, factor, max_terms, init_value);
0060 }
0061
0062 template <class Functor, class U>
0063 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, int bits, boost::math::uintmax_t& max_terms, const U& init_value) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0064 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0065 && noexcept(std::declval<Functor>()())
0066 #endif
0067 )
0068 {
0069    BOOST_MATH_STD_USING
0070    typedef typename Functor::result_type result_type;
0071    result_type factor = ldexp(result_type(1), 1 - bits);
0072    return sum_series(func, factor, max_terms, init_value);
0073 }
0074
0075 template <class Functor>
0076 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0077 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0078 && noexcept(std::declval<Functor>()())
0079 #endif
0080 )
0081 {
0082    BOOST_MATH_STD_USING
0083    typedef typename Functor::result_type result_type;
0084    boost::math::uintmax_t iters = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
0085    result_type init_val = 0;
0086    return sum_series(func, bits, iters, init_val);
0087 }
0088
0089 template <class Functor>
0090 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, int bits, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0091 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0092 && noexcept(std::declval<Functor>()())
0093 #endif
0094 )
0095 {
0096    BOOST_MATH_STD_USING
0097    typedef typename Functor::result_type result_type;
0098    result_type init_val = 0;
0099    return sum_series(func, bits, max_terms, init_val);
0100 }
0101
0102 template <class Functor, class U>
0103 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0104 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0105 && noexcept(std::declval<Functor>()())
0106 #endif
0107 )
0108 {
0109    BOOST_MATH_STD_USING
0110    boost::math::uintmax_t iters = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
0111    return sum_series(func, bits, iters, init_value);
0112 }
0113 //
0114 // Checked summation:
0115 //
0116 template <class Functor, class U, class V>
0117 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type checked_sum_series(Functor& func, const U& factor, boost::math::uintmax_t& max_terms, const V& init_value, V& norm) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0118 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0119 && noexcept(std::declval<Functor>()())
0120 #endif
0121 )
0122 {
0123    BOOST_MATH_STD_USING
0124
0125    typedef typename Functor::result_type result_type;
0126
0127    boost::math::uintmax_t counter = max_terms;
0128
0129    result_type result = init_value;
0130    result_type next_term;
0131    do {
0132       next_term = func();
0133       result += next_term;
0134       norm += fabs(next_term);
0135    } while ((abs(factor * result) < abs(next_term)) && --counter);
0136
0137    // set max_terms to the actual number of terms of the series evaluated:
0138    max_terms = max_terms - counter;
0139
0140    return result;
0141 }
0142
0143
0144 //
0145 // Algorithm kahan_sum_series invokes Functor func until the N'th
0146 // term is too small to have any effect on the total, the terms
0147 // are added using the Kahan summation method.
0148 //
0149 // CAUTION: Optimizing compilers combined with extended-precision
0150 // machine registers conspire to render this algorithm partly broken:
0151 // double rounding of intermediate terms (first to a long double machine
0152 // register, and then to a double result) cause the rounding error computed
0153 // by the algorithm to be off by up to 1ulp.  However this occurs rarely, and
0154 // in any case the result is still much better than a naive summation.
0155 //
0156 template <class Functor>
0157 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type kahan_sum_series(Functor& func, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0158 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0159 && noexcept(std::declval<Functor>()())
0160 #endif
0161 )
0162 {
0163    BOOST_MATH_STD_USING
0164
0165    typedef typename Functor::result_type result_type;
0166
0167    result_type factor = pow(result_type(2), result_type(bits));
0168    result_type result = func();
0169    result_type next_term, y, t;
0170    result_type carry = 0;
0171    do{
0172       next_term = func();
0173       y = next_term - carry;
0174       t = result + y;
0175       carry = t - result;
0176       carry -= y;
0177       result = t;
0178    }
0179    while(fabs(result) < fabs(factor * next_term));
0180    return result;
0181 }
0182
0183 template <class Functor>
0184 BOOST_MATH_GPU_ENABLED inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
0185 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0186 && noexcept(std::declval<Functor>()())
0187 #endif
0188 )
0189 {
0190    BOOST_MATH_STD_USING
0191
0192    typedef typename Functor::result_type result_type;
0193
0194    boost::math::uintmax_t counter = max_terms;
0195
0196    result_type factor = ldexp(result_type(1), bits);
0197    result_type result = func();
0198    result_type next_term, y, t;
0199    result_type carry = 0;
0200    do{
0201       next_term = func();
0202       y = next_term - carry;
0203       t = result + y;
0204       carry = t - result;
0205       carry -= y;
0206       result = t;
0207    }
0208    while((fabs(result) < fabs(factor * next_term)) && --counter);
0209
0210    // set max_terms to the actual number of terms of the series evaluated:
0211    max_terms = max_terms - counter;
0212
0213    return result;
0214 }
0215
0216 } // namespace tools
0217 } // namespace math
0218 } // namespace boost
0219
0220 #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED
0221