222 lines · plain
1// (C) Copyright John Maddock 2005-2006.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED7#define BOOST_MATH_TOOLS_SERIES_INCLUDED8 9#ifdef _MSC_VER10#pragma once11#endif12 13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/numeric_limits.hpp>16#include <boost/math/tools/cstdint.hpp>17#include <boost/math/tools/type_traits.hpp>18 19namespace boost{ namespace math{ namespace tools{20 21//22// Simple series summation come first:23//24template <class Functor, class U, class V>25BOOST_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)26#ifndef BOOST_MATH_HAS_GPU_SUPPORT27&& noexcept(std::declval<Functor>()())28#endif29)30{31 BOOST_MATH_STD_USING32 33 typedef typename Functor::result_type result_type;34 35 boost::math::uintmax_t counter = max_terms;36 37 result_type result = init_value;38 result_type next_term;39 do{40 next_term = func();41 result += next_term;42 }43 while((abs(factor * result) < abs(next_term)) && --counter);44 45 // set max_terms to the actual number of terms of the series evaluated:46 max_terms = max_terms - counter;47 48 return result;49}50 51template <class Functor, class U>52BOOST_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)53#ifndef BOOST_MATH_HAS_GPU_SUPPORT54&& noexcept(std::declval<Functor>()())55#endif56)57{58 typename Functor::result_type init_value = 0;59 return sum_series(func, factor, max_terms, init_value);60}61 62template <class Functor, class U>63BOOST_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)64#ifndef BOOST_MATH_HAS_GPU_SUPPORT65&& noexcept(std::declval<Functor>()())66#endif67)68{69 BOOST_MATH_STD_USING70 typedef typename Functor::result_type result_type;71 result_type factor = ldexp(result_type(1), 1 - bits);72 return sum_series(func, factor, max_terms, init_value);73}74 75template <class Functor>76BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type) 77#ifndef BOOST_MATH_HAS_GPU_SUPPORT78&& noexcept(std::declval<Functor>()())79#endif80)81{82 BOOST_MATH_STD_USING83 typedef typename Functor::result_type result_type;84 boost::math::uintmax_t iters = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();85 result_type init_val = 0;86 return sum_series(func, bits, iters, init_val);87}88 89template <class Functor>90BOOST_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)91#ifndef BOOST_MATH_HAS_GPU_SUPPORT92&& noexcept(std::declval<Functor>()())93#endif94)95{96 BOOST_MATH_STD_USING97 typedef typename Functor::result_type result_type;98 result_type init_val = 0;99 return sum_series(func, bits, max_terms, init_val);100}101 102template <class Functor, class U>103BOOST_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)104#ifndef BOOST_MATH_HAS_GPU_SUPPORT105&& noexcept(std::declval<Functor>()())106#endif107)108{109 BOOST_MATH_STD_USING110 boost::math::uintmax_t iters = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();111 return sum_series(func, bits, iters, init_value);112}113//114// Checked summation:115//116template <class Functor, class U, class V>117BOOST_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)118#ifndef BOOST_MATH_HAS_GPU_SUPPORT119&& noexcept(std::declval<Functor>()())120#endif121)122{123 BOOST_MATH_STD_USING124 125 typedef typename Functor::result_type result_type;126 127 boost::math::uintmax_t counter = max_terms;128 129 result_type result = init_value;130 result_type next_term;131 do {132 next_term = func();133 result += next_term;134 norm += fabs(next_term);135 } while ((abs(factor * result) < abs(next_term)) && --counter);136 137 // set max_terms to the actual number of terms of the series evaluated:138 max_terms = max_terms - counter;139 140 return result;141}142 143 144//145// Algorithm kahan_sum_series invokes Functor func until the N'th146// term is too small to have any effect on the total, the terms147// are added using the Kahan summation method.148//149// CAUTION: Optimizing compilers combined with extended-precision150// machine registers conspire to render this algorithm partly broken:151// double rounding of intermediate terms (first to a long double machine152// register, and then to a double result) cause the rounding error computed153// by the algorithm to be off by up to 1ulp. However this occurs rarely, and154// in any case the result is still much better than a naive summation.155//156template <class Functor>157BOOST_MATH_GPU_ENABLED inline typename Functor::result_type kahan_sum_series(Functor& func, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)158#ifndef BOOST_MATH_HAS_GPU_SUPPORT159&& noexcept(std::declval<Functor>()())160#endif161)162{163 BOOST_MATH_STD_USING164 165 typedef typename Functor::result_type result_type;166 167 result_type factor = pow(result_type(2), result_type(bits));168 result_type result = func();169 result_type next_term, y, t;170 result_type carry = 0;171 do{172 next_term = func();173 y = next_term - carry;174 t = result + y;175 carry = t - result;176 carry -= y;177 result = t;178 }179 while(fabs(result) < fabs(factor * next_term));180 return result;181}182 183template <class Functor>184BOOST_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)185#ifndef BOOST_MATH_HAS_GPU_SUPPORT186&& noexcept(std::declval<Functor>()())187#endif188)189{190 BOOST_MATH_STD_USING191 192 typedef typename Functor::result_type result_type;193 194 boost::math::uintmax_t counter = max_terms;195 196 result_type factor = ldexp(result_type(1), bits);197 result_type result = func();198 result_type next_term, y, t;199 result_type carry = 0;200 do{201 next_term = func();202 y = next_term - carry;203 t = result + y;204 carry = t - result;205 carry -= y;206 result = t;207 }208 while((fabs(result) < fabs(factor * next_term)) && --counter);209 210 // set max_terms to the actual number of terms of the series evaluated:211 max_terms = max_terms - counter;212 213 return result;214}215 216} // namespace tools217} // namespace math218} // namespace boost219 220#endif // BOOST_MATH_TOOLS_SERIES_INCLUDED221 222