math/special_functions/bessel_iterators.hpp

0001
0002 ///////////////////////////////////////////////////////////////////////////////
0003 //  Copyright 2018 John Maddock
0004 //  Distributed under the Boost
0005 //  Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_MATH_BESSEL_ITERATORS_HPP
0009 #define BOOST_MATH_BESSEL_ITERATORS_HPP
0010
0011 #include <boost/math/tools/recurrence.hpp>
0012 #include <boost/math/special_functions/bessel.hpp>
0013
0014 namespace boost {
0015    namespace math {
0016       namespace detail {
0017
0018          template <class T>
0019          struct bessel_jy_recurrence
0020          {
0021             bessel_jy_recurrence(T v, T z) : v(v), z(z) {}
0022             boost::math::tuple<T, T, T> operator()(int k)
0023             {
0024                return boost::math::tuple<T, T, T>(1, -2 * (v + k) / z, 1);
0025             }
0026
0027             T v, z;
0028          };
0029          template <class T>
0030          struct bessel_ik_recurrence
0031          {
0032             bessel_ik_recurrence(T v, T z) : v(v), z(z) {}
0033             boost::math::tuple<T, T, T> operator()(int k)
0034             {
0035                return boost::math::tuple<T, T, T>(1, -2 * (v + k) / z, -1);
0036             }
0037
0038             T v, z;
0039          };
0040       } // namespace detail
0041
0042       template <class T, class Policy = boost::math::policies::policy<> >
0043       struct bessel_j_backwards_iterator
0044       {
0045          typedef std::ptrdiff_t difference_type;
0046          typedef T value_type;
0047          typedef T* pointer;
0048          typedef T& reference;
0049          typedef std::input_iterator_tag iterator_category;
0050
0051          bessel_j_backwards_iterator(const T& v, const T& x)
0052             : it(detail::bessel_jy_recurrence<T>(v, x), boost::math::cyl_bessel_j(v, x, Policy()))
0053          {
0054             if(v < 0)
0055                boost::math::policies::raise_domain_error("bessel_j_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, Policy());
0056          }
0057
0058          bessel_j_backwards_iterator(const T& v, const T& x, const T& J_v)
0059             : it(detail::bessel_jy_recurrence<T>(v, x), J_v)
0060          {
0061             if(v < 0)
0062                boost::math::policies::raise_domain_error("bessel_j_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, Policy());
0063          }
0064          bessel_j_backwards_iterator(const T& v, const T& x, const T& J_v_plus_1, const T& J_v)
0065             : it(detail::bessel_jy_recurrence<T>(v, x), J_v_plus_1, J_v)
0066          {
0067             if (v < -1)
0068                boost::math::policies::raise_domain_error("bessel_j_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, Policy());
0069          }
0070
0071          bessel_j_backwards_iterator& operator++()
0072          {
0073             ++it;
0074             return *this;
0075          }
0076
0077          bessel_j_backwards_iterator operator++(int)
0078          {
0079             bessel_j_backwards_iterator t(*this);
0080             ++(*this);
0081             return t;
0082          }
0083
0084          T operator*() { return *it; }
0085
0086       private:
0087          boost::math::tools::backward_recurrence_iterator< detail::bessel_jy_recurrence<T> > it;
0088       };
0089
0090       template <class T, class Policy = boost::math::policies::policy<> >
0091       struct bessel_i_backwards_iterator
0092       {
0093          typedef std::ptrdiff_t difference_type;
0094          typedef T value_type;
0095          typedef T* pointer;
0096          typedef T& reference;
0097          typedef std::input_iterator_tag iterator_category;
0098
0099          bessel_i_backwards_iterator(const T& v, const T& x)
0100             : it(detail::bessel_ik_recurrence<T>(v, x), boost::math::cyl_bessel_i(v, x, Policy()))
0101          {
0102             if(v < -1)
0103                boost::math::policies::raise_domain_error("bessel_i_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, Policy());
0104          }
0105          bessel_i_backwards_iterator(const T& v, const T& x, const T& I_v)
0106             : it(detail::bessel_ik_recurrence<T>(v, x), I_v)
0107          {
0108             if(v < -1)
0109                boost::math::policies::raise_domain_error("bessel_i_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, Policy());
0110          }
0111          bessel_i_backwards_iterator(const T& v, const T& x, const T& I_v_plus_1, const T& I_v)
0112             : it(detail::bessel_ik_recurrence<T>(v, x), I_v_plus_1, I_v)
0113          {
0114             if(v < -1)
0115                boost::math::policies::raise_domain_error("bessel_i_backwards_iterator<%1%>", "Order must be > 0 stable backwards recurrence but got %1%", v, Policy());
0116          }
0117
0118          bessel_i_backwards_iterator& operator++()
0119          {
0120             ++it;
0121             return *this;
0122          }
0123
0124          bessel_i_backwards_iterator operator++(int)
0125          {
0126             bessel_i_backwards_iterator t(*this);
0127             ++(*this);
0128             return t;
0129          }
0130
0131          T operator*() { return *it; }
0132
0133       private:
0134          boost::math::tools::backward_recurrence_iterator< detail::bessel_ik_recurrence<T> > it;
0135       };
0136
0137       template <class T, class Policy = boost::math::policies::policy<> >
0138       struct bessel_i_forwards_iterator
0139       {
0140          typedef std::ptrdiff_t difference_type;
0141          typedef T value_type;
0142          typedef T* pointer;
0143          typedef T& reference;
0144          typedef std::input_iterator_tag iterator_category;
0145
0146          bessel_i_forwards_iterator(const T& v, const T& x)
0147             : it(detail::bessel_ik_recurrence<T>(v, x), boost::math::cyl_bessel_i(v, x, Policy()))
0148          {
0149             if(v > 1)
0150                boost::math::policies::raise_domain_error("bessel_i_forwards_iterator<%1%>", "Order must be < 0 stable forwards recurrence but got %1%", v, Policy());
0151          }
0152          bessel_i_forwards_iterator(const T& v, const T& x, const T& I_v)
0153             : it(detail::bessel_ik_recurrence<T>(v, x), I_v)
0154          {
0155             if (v > 1)
0156                boost::math::policies::raise_domain_error("bessel_i_forwards_iterator<%1%>", "Order must be < 0 stable forwards recurrence but got %1%", v, Policy());
0157          }
0158          bessel_i_forwards_iterator(const T& v, const T& x, const T& I_v_minus_1, const T& I_v)
0159             : it(detail::bessel_ik_recurrence<T>(v, x), I_v_minus_1, I_v)
0160          {
0161             if (v > 1)
0162                boost::math::policies::raise_domain_error("bessel_i_forwards_iterator<%1%>", "Order must be < 0 stable forwards recurrence but got %1%", v, Policy());
0163          }
0164
0165          bessel_i_forwards_iterator& operator++()
0166          {
0167             ++it;
0168             return *this;
0169          }
0170
0171          bessel_i_forwards_iterator operator++(int)
0172          {
0173             bessel_i_forwards_iterator t(*this);
0174             ++(*this);
0175             return t;
0176          }
0177
0178          T operator*() { return *it; }
0179
0180       private:
0181          boost::math::tools::forward_recurrence_iterator< detail::bessel_ik_recurrence<T> > it;
0182       };
0183
0184    }
0185 } // namespaces
0186
0187 #endif // BOOST_MATH_BESSEL_ITERATORS_HPP