series.hpp 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158
  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. #ifndef BOOST_MATH_TOOLS_SERIES_INCLUDED
  6. #define BOOST_MATH_TOOLS_SERIES_INCLUDED
  7. #ifdef _MSC_VER
  8. #pragma once
  9. #endif
  10. #include <boost/config/no_tr1/cmath.hpp>
  11. #include <boost/cstdint.hpp>
  12. #include <boost/limits.hpp>
  13. #include <boost/math/tools/config.hpp>
  14. namespace boost{ namespace math{ namespace tools{
  15. //
  16. // Simple series summation come first:
  17. //
  18. template <class Functor, class U, class V>
  19. inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms, const V& init_value) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  20. {
  21. BOOST_MATH_STD_USING
  22. typedef typename Functor::result_type result_type;
  23. boost::uintmax_t counter = max_terms;
  24. result_type result = init_value;
  25. result_type next_term;
  26. do{
  27. next_term = func();
  28. result += next_term;
  29. }
  30. while((fabs(factor * result) < fabs(next_term)) && --counter);
  31. // set max_terms to the actual number of terms of the series evaluated:
  32. max_terms = max_terms - counter;
  33. return result;
  34. }
  35. template <class Functor, class U>
  36. inline typename Functor::result_type sum_series(Functor& func, const U& factor, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  37. {
  38. typename Functor::result_type init_value = 0;
  39. return sum_series(func, factor, max_terms, init_value);
  40. }
  41. template <class Functor, class U>
  42. inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  43. {
  44. BOOST_MATH_STD_USING
  45. typedef typename Functor::result_type result_type;
  46. result_type factor = ldexp(result_type(1), 1 - bits);
  47. return sum_series(func, factor, max_terms, init_value);
  48. }
  49. template <class Functor>
  50. inline typename Functor::result_type sum_series(Functor& func, int bits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  51. {
  52. BOOST_MATH_STD_USING
  53. typedef typename Functor::result_type result_type;
  54. boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
  55. result_type init_val = 0;
  56. return sum_series(func, bits, iters, init_val);
  57. }
  58. template <class Functor>
  59. inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  60. {
  61. BOOST_MATH_STD_USING
  62. typedef typename Functor::result_type result_type;
  63. result_type init_val = 0;
  64. return sum_series(func, bits, max_terms, init_val);
  65. }
  66. template <class Functor, class U>
  67. inline typename Functor::result_type sum_series(Functor& func, int bits, const U& init_value) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  68. {
  69. BOOST_MATH_STD_USING
  70. boost::uintmax_t iters = (std::numeric_limits<boost::uintmax_t>::max)();
  71. return sum_series(func, bits, iters, init_value);
  72. }
  73. //
  74. // Algorithm kahan_sum_series invokes Functor func until the N'th
  75. // term is too small to have any effect on the total, the terms
  76. // are added using the Kahan summation method.
  77. //
  78. // CAUTION: Optimizing compilers combined with extended-precision
  79. // machine registers conspire to render this algorithm partly broken:
  80. // double rounding of intermediate terms (first to a long double machine
  81. // register, and then to a double result) cause the rounding error computed
  82. // by the algorithm to be off by up to 1ulp. However this occurs rarely, and
  83. // in any case the result is still much better than a naive summation.
  84. //
  85. template <class Functor>
  86. inline typename Functor::result_type kahan_sum_series(Functor& func, int bits) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  87. {
  88. BOOST_MATH_STD_USING
  89. typedef typename Functor::result_type result_type;
  90. result_type factor = pow(result_type(2), bits);
  91. result_type result = func();
  92. result_type next_term, y, t;
  93. result_type carry = 0;
  94. do{
  95. next_term = func();
  96. y = next_term - carry;
  97. t = result + y;
  98. carry = t - result;
  99. carry -= y;
  100. result = t;
  101. }
  102. while(fabs(result) < fabs(factor * next_term));
  103. return result;
  104. }
  105. template <class Functor>
  106. inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(typename Functor::result_type) && noexcept(std::declval<Functor>()()))
  107. {
  108. BOOST_MATH_STD_USING
  109. typedef typename Functor::result_type result_type;
  110. boost::uintmax_t counter = max_terms;
  111. result_type factor = ldexp(result_type(1), bits);
  112. result_type result = func();
  113. result_type next_term, y, t;
  114. result_type carry = 0;
  115. do{
  116. next_term = func();
  117. y = next_term - carry;
  118. t = result + y;
  119. carry = t - result;
  120. carry -= y;
  121. result = t;
  122. }
  123. while((fabs(result) < fabs(factor * next_term)) && --counter);
  124. // set max_terms to the actual number of terms of the series evaluated:
  125. max_terms = max_terms - counter;
  126. return result;
  127. }
  128. } // namespace tools
  129. } // namespace math
  130. } // namespace boost
  131. #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED