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