series.hpp 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222
  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/math/tools/config.hpp>
  11. #include <boost/math/tools/numeric_limits.hpp>
  12. #include <boost/math/tools/cstdint.hpp>
  13. #include <boost/math/tools/type_traits.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. 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)
  20. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  21. && noexcept(std::declval<Functor>()())
  22. #endif
  23. )
  24. {
  25. BOOST_MATH_STD_USING
  26. typedef typename Functor::result_type result_type;
  27. boost::math::uintmax_t counter = max_terms;
  28. result_type result = init_value;
  29. result_type next_term;
  30. do{
  31. next_term = func();
  32. result += next_term;
  33. }
  34. while((abs(factor * result) < abs(next_term)) && --counter);
  35. // set max_terms to the actual number of terms of the series evaluated:
  36. max_terms = max_terms - counter;
  37. return result;
  38. }
  39. template <class Functor, class U>
  40. 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)
  41. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  42. && noexcept(std::declval<Functor>()())
  43. #endif
  44. )
  45. {
  46. typename Functor::result_type init_value = 0;
  47. return sum_series(func, factor, max_terms, init_value);
  48. }
  49. template <class Functor, class U>
  50. 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)
  51. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  52. && noexcept(std::declval<Functor>()())
  53. #endif
  54. )
  55. {
  56. BOOST_MATH_STD_USING
  57. typedef typename Functor::result_type result_type;
  58. result_type factor = ldexp(result_type(1), 1 - bits);
  59. return sum_series(func, factor, max_terms, init_value);
  60. }
  61. template <class Functor>
  62. BOOST_MATH_GPU_ENABLED inline typename Functor::result_type sum_series(Functor& func, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename Functor::result_type)
  63. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  64. && noexcept(std::declval<Functor>()())
  65. #endif
  66. )
  67. {
  68. BOOST_MATH_STD_USING
  69. typedef typename Functor::result_type result_type;
  70. boost::math::uintmax_t iters = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
  71. result_type init_val = 0;
  72. return sum_series(func, bits, iters, init_val);
  73. }
  74. template <class Functor>
  75. 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)
  76. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  77. && noexcept(std::declval<Functor>()())
  78. #endif
  79. )
  80. {
  81. BOOST_MATH_STD_USING
  82. typedef typename Functor::result_type result_type;
  83. result_type init_val = 0;
  84. return sum_series(func, bits, max_terms, init_val);
  85. }
  86. template <class Functor, class U>
  87. 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)
  88. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  89. && noexcept(std::declval<Functor>()())
  90. #endif
  91. )
  92. {
  93. BOOST_MATH_STD_USING
  94. boost::math::uintmax_t iters = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
  95. return sum_series(func, bits, iters, init_value);
  96. }
  97. //
  98. // Checked summation:
  99. //
  100. template <class Functor, class U, class V>
  101. 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)
  102. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  103. && noexcept(std::declval<Functor>()())
  104. #endif
  105. )
  106. {
  107. BOOST_MATH_STD_USING
  108. typedef typename Functor::result_type result_type;
  109. boost::math::uintmax_t counter = max_terms;
  110. result_type result = init_value;
  111. result_type next_term;
  112. do {
  113. next_term = func();
  114. result += next_term;
  115. norm += fabs(next_term);
  116. } while ((abs(factor * result) < abs(next_term)) && --counter);
  117. // set max_terms to the actual number of terms of the series evaluated:
  118. max_terms = max_terms - counter;
  119. return result;
  120. }
  121. //
  122. // Algorithm kahan_sum_series invokes Functor func until the N'th
  123. // term is too small to have any effect on the total, the terms
  124. // are added using the Kahan summation method.
  125. //
  126. // CAUTION: Optimizing compilers combined with extended-precision
  127. // machine registers conspire to render this algorithm partly broken:
  128. // double rounding of intermediate terms (first to a long double machine
  129. // register, and then to a double result) cause the rounding error computed
  130. // by the algorithm to be off by up to 1ulp. However this occurs rarely, and
  131. // in any case the result is still much better than a naive summation.
  132. //
  133. template <class Functor>
  134. 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)
  135. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  136. && noexcept(std::declval<Functor>()())
  137. #endif
  138. )
  139. {
  140. BOOST_MATH_STD_USING
  141. typedef typename Functor::result_type result_type;
  142. result_type factor = pow(result_type(2), result_type(bits));
  143. result_type result = func();
  144. result_type next_term, y, t;
  145. result_type carry = 0;
  146. do{
  147. next_term = func();
  148. y = next_term - carry;
  149. t = result + y;
  150. carry = t - result;
  151. carry -= y;
  152. result = t;
  153. }
  154. while(fabs(result) < fabs(factor * next_term));
  155. return result;
  156. }
  157. template <class Functor>
  158. 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)
  159. #ifndef BOOST_MATH_HAS_GPU_SUPPORT
  160. && noexcept(std::declval<Functor>()())
  161. #endif
  162. )
  163. {
  164. BOOST_MATH_STD_USING
  165. typedef typename Functor::result_type result_type;
  166. boost::math::uintmax_t counter = max_terms;
  167. result_type factor = ldexp(result_type(1), 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)) && --counter);
  180. // set max_terms to the actual number of terms of the series evaluated:
  181. max_terms = max_terms - counter;
  182. return result;
  183. }
  184. } // namespace tools
  185. } // namespace math
  186. } // namespace boost
  187. #endif // BOOST_MATH_TOOLS_SERIES_INCLUDED