annotate DEPENDENCIES/generic/include/boost/math/tools/series.hpp @ 133:4acb5d8d80b6 tip

Don't fail environmental check if README.md exists (but .txt and no-suffix don't)
author Chris Cannam
date Tue, 30 Jul 2019 12:25:44 +0100
parents 2665513ce2d3
children
rev   line source
Chris@16 1 // (C) Copyright John Maddock 2005-2006.
Chris@16 2 // Use, modification and distribution are subject to the
Chris@16 3 // Boost Software License, Version 1.0. (See accompanying file
Chris@16 4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16 5
Chris@16 6 #ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED
Chris@16 7 #define BOOST_MATH_TOOLS_SERIES_INCLUDED
Chris@16 8
Chris@16 9 #ifdef _MSC_VER
Chris@16 10 #pragma once
Chris@16 11 #endif
Chris@16 12
Chris@16 13 #include <boost/config/no_tr1/cmath.hpp>
Chris@16 14 #include <boost/cstdint.hpp>
Chris@16 15 #include <boost/limits.hpp>
Chris@16 16 #include <boost/math/tools/config.hpp>
Chris@16 17
Chris@16 18 namespace boost{ namespace math{ namespace tools{
Chris@16 19
Chris@16 20 //
Chris@16 21 // Simple series summation come first:
Chris@16 22 //
Chris@16 23 template <class Functor, class U, class V>
Chris@16 24 inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value)
Chris@16 25 {
Chris@16 26 BOOST_MATH_STD_USING
Chris@16 27
Chris@16 28 typedef typename Functor::result_type result_type;
Chris@16 29
Chris@16 30 boost::uintmax_t counter = max_terms;
Chris@16 31
Chris@16 32 result_type result = init_value;
Chris@16 33 result_type next_term;
Chris@16 34 do{
Chris@16 35 next_term = func();
Chris@16 36 result += next_term;
Chris@16 37 }
Chris@16 38 while((fabs(factor * result) < fabs(next_term)) && --counter);
Chris@16 39
Chris@16 40 // set max_terms to the actual number of terms of the series evaluated:
Chris@16 41 max_terms = max_terms - counter;
Chris@16 42
Chris@16 43 return result;
Chris@16 44 }
Chris@16 45
Chris@16 46 template <class Functor, class U>
Chris@16 47 inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms)
Chris@16 48 {
Chris@16 49 typename Functor::result_type init_value = 0;
Chris@16 50 return sum_series(func, factor, max_terms, init_value);
Chris@16 51 }
Chris@16 52
Chris@16 53 template <class Functor, class U>
Chris@16 54 inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value)
Chris@16 55 {
Chris@16 56 BOOST_MATH_STD_USING
Chris@16 57 typedef typename Functor::result_type result_type;
Chris@16 58 result_type factor = ldexp(result_type(1), 1 - bits);
Chris@16 59 return sum_series(func, factor, max_terms, init_value);
Chris@16 60 }
Chris@16 61
Chris@16 62 template <class Functor>
Chris@16 63 inline typename Functor::result_type sum_series(Functor& func, int bits)
Chris@16 64 {
Chris@16 65 BOOST_MATH_STD_USING
Chris@16 66 typedef typename Functor::result_type result_type;
Chris@16 67 boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16 68 result_type init_val = 0;
Chris@16 69 return sum_series(func, bits, iters, init_val);
Chris@16 70 }
Chris@16 71
Chris@16 72 template <class Functor>
Chris@16 73 inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
Chris@16 74 {
Chris@16 75 BOOST_MATH_STD_USING
Chris@16 76 typedef typename Functor::result_type result_type;
Chris@16 77 result_type init_val = 0;
Chris@16 78 return sum_series(func, bits, max_terms, init_val);
Chris@16 79 }
Chris@16 80
Chris@16 81 template <class Functor, class U>
Chris@16 82 inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value)
Chris@16 83 {
Chris@16 84 BOOST_MATH_STD_USING
Chris@16 85 boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
Chris@16 86 return sum_series(func, bits, iters, init_value);
Chris@16 87 }
Chris@16 88
Chris@16 89 //
Chris@16 90 // Algorithm kahan_sum_series invokes Functor func until the N'th
Chris@16 91 // term is too small to have any effect on the total, the terms
Chris@16 92 // are added using the Kahan summation method.
Chris@16 93 //
Chris@16 94 // CAUTION: Optimizing compilers combined with extended-precision
Chris@16 95 // machine registers conspire to render this algorithm partly broken:
Chris@16 96 // double rounding of intermediate terms (first to a long double machine
Chris@16 97 // register, and then to a double result) cause the rounding error computed
Chris@16 98 // by the algorithm to be off by up to 1ulp. However this occurs rarely, and
Chris@16 99 // in any case the result is still much better than a naive summation.
Chris@16 100 //
Chris@16 101 template <class Functor>
Chris@16 102 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits)
Chris@16 103 {
Chris@16 104 BOOST_MATH_STD_USING
Chris@16 105
Chris@16 106 typedef typename Functor::result_type result_type;
Chris@16 107
Chris@16 108 result_type factor = pow(result_type(2), bits);
Chris@16 109 result_type result = func();
Chris@16 110 result_type next_term, y, t;
Chris@16 111 result_type carry = 0;
Chris@16 112 do{
Chris@16 113 next_term = func();
Chris@16 114 y = next_term - carry;
Chris@16 115 t = result + y;
Chris@16 116 carry = t - result;
Chris@16 117 carry -= y;
Chris@16 118 result = t;
Chris@16 119 }
Chris@16 120 while(fabs(result) < fabs(factor * next_term));
Chris@16 121 return result;
Chris@16 122 }
Chris@16 123
Chris@16 124 template <class Functor>
Chris@16 125 inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
Chris@16 126 {
Chris@16 127 BOOST_MATH_STD_USING
Chris@16 128
Chris@16 129 typedef typename Functor::result_type result_type;
Chris@16 130
Chris@16 131 boost::uintmax_t counter = max_terms;
Chris@16 132
Chris@16 133 result_type factor = ldexp(result_type(1), bits);
Chris@16 134 result_type result = func();
Chris@16 135 result_type next_term, y, t;
Chris@16 136 result_type carry = 0;
Chris@16 137 do{
Chris@16 138 next_term = func();
Chris@16 139 y = next_term - carry;
Chris@16 140 t = result + y;
Chris@16 141 carry = t - result;
Chris@16 142 carry -= y;
Chris@16 143 result = t;
Chris@16 144 }
Chris@16 145 while((fabs(result) < fabs(factor * next_term)) && --counter);
Chris@16 146
Chris@16 147 // set max_terms to the actual number of terms of the series evaluated:
Chris@16 148 max_terms = max_terms - counter;
Chris@16 149
Chris@16 150 return result;
Chris@16 151 }
Chris@16 152
Chris@16 153 } // namespace tools
Chris@16 154 } // namespace math
Chris@16 155 } // namespace boost
Chris@16 156
Chris@16 157 #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED
Chris@16 158