bessel_iterators.hpp 6.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // Copyright 2018 John Maddock
  3. // Distributed under the Boost
  4. // Software License, Version 1.0. (See accompanying file
  5. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. #ifndef BOOST_MATH_BESSEL_ITERATORS_HPP
  7. #define BOOST_MATH_BESSEL_ITERATORS_HPP
  8. #include <boost/math/tools/recurrence.hpp>
  9. namespace boost {
  10. namespace math {
  11. namespace detail {
  12. template <class T>
  13. struct bessel_jy_recurrence
  14. {
  15. bessel_jy_recurrence(T v, T z) : v(v), z(z) {}
  16. boost::math::tuple<T, T, T> operator()(int k)
  17. {
  18. return boost::math::tuple<T, T, T>(1, -2 * (v + k) / z, 1);
  19. }
  20. T v, z;
  21. };
  22. template <class T>
  23. struct bessel_ik_recurrence
  24. {
  25. bessel_ik_recurrence(T v, T z) : v(v), z(z) {}
  26. boost::math::tuple<T, T, T> operator()(int k)
  27. {
  28. return boost::math::tuple<T, T, T>(1, -2 * (v + k) / z, -1);
  29. }
  30. T v, z;
  31. };
  32. } // namespace detail
  33. template <class T>
  34. struct bessel_j_backwards_iterator
  35. {
  36. typedef std::ptrdiff_t difference_type;
  37. typedef T value_type;
  38. typedef T* pointer;
  39. typedef T& reference;
  40. typedef std::input_iterator_tag iterator_category;
  41. bessel_j_backwards_iterator(const T& v, const T& x)
  42. : it(detail::bessel_jy_recurrence<T>(v, x), boost::math::cyl_bessel_j(v, x))
  43. {
  44. if(v < 0)
  45. boost::math::policies::raise_domain_error("bessel_j_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, boost::math::policies::policy<>());
  46. }
  47. bessel_j_backwards_iterator(const T& v, const T& x, const T& J_v)
  48. : it(detail::bessel_jy_recurrence<T>(v, x), J_v)
  49. {
  50. if(v < 0)
  51. boost::math::policies::raise_domain_error("bessel_j_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, boost::math::policies::policy<>());
  52. }
  53. bessel_j_backwards_iterator(const T& v, const T& x, const T& J_v_plus_1, const T& J_v)
  54. : it(detail::bessel_jy_recurrence<T>(v, x), J_v_plus_1, J_v)
  55. {
  56. if (v < -1)
  57. boost::math::policies::raise_domain_error("bessel_j_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, boost::math::policies::policy<>());
  58. }
  59. bessel_j_backwards_iterator& operator++()
  60. {
  61. ++it;
  62. return *this;
  63. }
  64. bessel_j_backwards_iterator operator++(int)
  65. {
  66. bessel_j_backwards_iterator t(*this);
  67. ++(*this);
  68. return t;
  69. }
  70. T operator*() { return *it; }
  71. private:
  72. boost::math::tools::backward_recurrence_iterator< detail::bessel_jy_recurrence<T> > it;
  73. };
  74. template <class T>
  75. struct bessel_i_backwards_iterator
  76. {
  77. typedef std::ptrdiff_t difference_type;
  78. typedef T value_type;
  79. typedef T* pointer;
  80. typedef T& reference;
  81. typedef std::input_iterator_tag iterator_category;
  82. bessel_i_backwards_iterator(const T& v, const T& x)
  83. : it(detail::bessel_ik_recurrence<T>(v, x), boost::math::cyl_bessel_i(v, x))
  84. {
  85. if(v < -1)
  86. boost::math::policies::raise_domain_error("bessel_i_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, boost::math::policies::policy<>());
  87. }
  88. bessel_i_backwards_iterator(const T& v, const T& x, const T& I_v)
  89. : it(detail::bessel_ik_recurrence<T>(v, x), I_v)
  90. {
  91. if(v < -1)
  92. boost::math::policies::raise_domain_error("bessel_i_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, boost::math::policies::policy<>());
  93. }
  94. bessel_i_backwards_iterator(const T& v, const T& x, const T& I_v_plus_1, const T& I_v)
  95. : it(detail::bessel_ik_recurrence<T>(v, x), I_v_plus_1, I_v)
  96. {
  97. if(v < -1)
  98. boost::math::policies::raise_domain_error("bessel_i_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, boost::math::policies::policy<>());
  99. }
  100. bessel_i_backwards_iterator& operator++()
  101. {
  102. ++it;
  103. return *this;
  104. }
  105. bessel_i_backwards_iterator operator++(int)
  106. {
  107. bessel_i_backwards_iterator t(*this);
  108. ++(*this);
  109. return t;
  110. }
  111. T operator*() { return *it; }
  112. private:
  113. boost::math::tools::backward_recurrence_iterator< detail::bessel_ik_recurrence<T> > it;
  114. };
  115. template <class T>
  116. struct bessel_i_forwards_iterator
  117. {
  118. typedef std::ptrdiff_t difference_type;
  119. typedef T value_type;
  120. typedef T* pointer;
  121. typedef T& reference;
  122. typedef std::input_iterator_tag iterator_category;
  123. bessel_i_forwards_iterator(const T& v, const T& x)
  124. : it(detail::bessel_ik_recurrence<T>(v, x), boost::math::cyl_bessel_i(v, x))
  125. {
  126. if(v > 1)
  127. boost::math::policies::raise_domain_error("bessel_i_forwards_iterator<%1%>", "Order must be < 0 stable forwards recurrence but got %1%", v, boost::math::policies::policy<>());
  128. }
  129. bessel_i_forwards_iterator(const T& v, const T& x, const T& I_v)
  130. : it(detail::bessel_ik_recurrence<T>(v, x), I_v)
  131. {
  132. if (v > 1)
  133. boost::math::policies::raise_domain_error("bessel_i_forwards_iterator<%1%>", "Order must be < 0 stable forwards recurrence but got %1%", v, boost::math::policies::policy<>());
  134. }
  135. bessel_i_forwards_iterator(const T& v, const T& x, const T& I_v_plus_1, const T& I_v)
  136. : it(detail::bessel_ik_recurrence<T>(v, x), I_v_plus_1, I_v)
  137. {
  138. if (v > 1)
  139. boost::math::policies::raise_domain_error("bessel_i_forwards_iterator<%1%>", "Order must be < 0 stable forwards recurrence but got %1%", v, boost::math::policies::policy<>());
  140. }
  141. bessel_i_forwards_iterator& operator++()
  142. {
  143. ++it;
  144. return *this;
  145. }
  146. bessel_i_forwards_iterator operator++(int)
  147. {
  148. bessel_i_forwards_iterator t(*this);
  149. ++(*this);
  150. return t;
  151. }
  152. T operator*() { return *it; }
  153. private:
  154. boost::math::tools::forward_recurrence_iterator< detail::bessel_ik_recurrence<T> > it;
  155. };
  156. }
  157. } // namespaces
  158. #endif // BOOST_MATH_BESSEL_ITERATORS_HPP