boost/math/special_functions/bessel_iterators.hpp

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