boost/math/special_functions/detail/hypergeometric_cf.hpp

   1
   2 ///////////////////////////////////////////////////////////////////////////////
   3 //  Copyright 2014 Anton Bikineev
   4 //  Copyright 2014 Christopher Kormanyos
   5 //  Copyright 2014 John Maddock
   6 //  Copyright 2014 Paul Bristow
   7 //  Distributed under the Boost
   8 //  Software License, Version 1.0. (See accompanying file
   9 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  10 //
  11 #ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP
  12 #define BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP
  13
  14   namespace boost { namespace math { namespace detail {
  15
  16   // primary template for term of continued fraction
  17   template <class T, unsigned p, unsigned q>
  18   struct hypergeometric_pFq_cf_term;
  19
  20   // partial specialization for 0F1
  21   template <class T>
  22   struct hypergeometric_pFq_cf_term<T, 0u, 1u>
  23   {
  24     typedef std::pair<T,T> result_type;
  25
  26     hypergeometric_pFq_cf_term(const T& b, const T& z):
  27       n(1), b(b), z(z),
  28       term(std::make_pair(T(0), T(1)))
  29     {
  30     }
  31
  32     result_type operator()()
  33     {
  34       const result_type result = term;
  35       ++b; ++n;
  36       numer = -(z / (b * n));
  37       term = std::make_pair(numer, 1 - numer);
  38       return result;
  39     }
  40
  41   private:
  42     unsigned n;
  43     T b;
  44     const T z;
  45     T numer;
  46     result_type term;
  47   };
  48
  49   // partial specialization for 1F0
  50   template <class T>
  51   struct hypergeometric_pFq_cf_term<T, 1u, 0u>
  52   {
  53     typedef std::pair<T,T> result_type;
  54
  55     hypergeometric_pFq_cf_term(const T& a, const T& z):
  56       n(1), a(a), z(z),
  57       term(std::make_pair(T(0), T(1)))
  58     {
  59     }
  60
  61     result_type operator()()
  62     {
  63       const result_type result = term;
  64       ++a; ++n;
  65       numer = -((a * z) / n);
  66       term = std::make_pair(numer, 1 - numer);
  67       return result;
  68     }
  69
  70   private:
  71     unsigned n;
  72     T a;
  73     const T z;
  74     T numer;
  75     result_type term;
  76   };
  77
  78   // partial specialization for 1F1
  79   template <class T>
  80   struct hypergeometric_pFq_cf_term<T, 1u, 1u>
  81   {
  82     typedef std::pair<T,T> result_type;
  83
  84     hypergeometric_pFq_cf_term(const T& a, const T& b, const T& z):
  85       n(1), a(a), b(b), z(z),
  86       term(std::make_pair(T(0), T(1)))
  87     {
  88     }
  89
  90     result_type operator()()
  91     {
  92       const result_type result = term;
  93       ++a; ++b; ++n;
  94       numer = -((a * z) / (b * n));
  95       term = std::make_pair(numer, 1 - numer);
  96       return result;
  97     }
  98
  99   private:
 100     unsigned n;
 101     T a, b;
 102     const T z;
 103     T numer;
 104     result_type term;
 105   };
 106
 107   // partial specialization for 1f2
 108   template <class T>
 109   struct hypergeometric_pFq_cf_term<T, 1u, 2u>
 110   {
 111     typedef std::pair<T,T> result_type;
 112
 113     hypergeometric_pFq_cf_term(const T& a, const T& b, const T& c, const T& z):
 114       n(1), a(a), b(b), c(c), z(z),
 115       term(std::make_pair(T(0), T(1)))
 116     {
 117     }
 118
 119     result_type operator()()
 120     {
 121       const result_type result = term;
 122       ++a; ++b; ++c; ++n;
 123       numer = -((a * z) / ((b * c) * n));
 124       term = std::make_pair(numer, 1 - numer);
 125       return result;
 126     }
 127
 128   private:
 129     unsigned n;
 130     T a, b, c;
 131     const T z;
 132     T numer;
 133     result_type term;
 134   };
 135
 136   // partial specialization for 2f1
 137   template <class T>
 138   struct hypergeometric_pFq_cf_term<T, 2u, 1u>
 139   {
 140     typedef std::pair<T,T> result_type;
 141
 142     hypergeometric_pFq_cf_term(const T& a, const T& b, const T& c, const T& z):
 143       n(1), a(a), b(b), c(c), z(z),
 144       term(std::make_pair(T(0), T(1)))
 145     {
 146     }
 147
 148     result_type operator()()
 149     {
 150       const result_type result = term;
 151       ++a; ++b; ++c; ++n;
 152       numer = -(((a * b) * z) / (c * n));
 153       term = std::make_pair(numer, 1 - numer);
 154       return result;
 155     }
 156
 157   private:
 158     unsigned n;
 159     T a, b, c;
 160     const T z;
 161     T numer;
 162     result_type term;
 163   };
 164
 165   template <class T, unsigned p, unsigned q, class Policy>
 166   inline T compute_cf_pFq(detail::hypergeometric_pFq_cf_term<T, p, q>& term, const Policy& pol)
 167   {
 168     BOOST_MATH_STD_USING
 169     boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 170     const T result = tools::continued_fraction_b(
 171       term,
 172       boost::math::policies::get_epsilon<T, Policy>(),
 173       max_iter);
 174     boost::math::policies::check_series_iterations<T>(
 175       "boost::math::hypergeometric_pFq_cf<%1%>(%1%,%1%,%1%)",
 176       max_iter,
 177       pol);
 178     return result;
 179   }
 180
 181   template <class T, class Policy>
 182   inline T hypergeometric_0F1_cf(const T& b, const T& z, const Policy& pol)
 183   {
 184     detail::hypergeometric_pFq_cf_term<T, 0u, 1u> f(b, z);
 185     T result = detail::compute_cf_pFq(f, pol);
 186     result = ((z / b) / result) + 1;
 187     return result;
 188   }
 189
 190   template <class T, class Policy>
 191   inline T hypergeometric_1F0_cf(const T& a, const T& z, const Policy& pol)
 192   {
 193     detail::hypergeometric_pFq_cf_term<T, 1u, 0u> f(a, z);
 194     T result = detail::compute_cf_pFq(f, pol);
 195     result = ((a * z) / result) + 1;
 196     return result;
 197   }
 198
 199   template <class T, class Policy>
 200   inline T hypergeometric_1F1_cf(const T& a, const T& b, const T& z, const Policy& pol)
 201   {
 202     detail::hypergeometric_pFq_cf_term<T, 1u, 1u> f(a, b, z);
 203     T result = detail::compute_cf_pFq(f, pol);
 204     result = (((a * z) / b) / result) + 1;
 205     return result;
 206   }
 207
 208   template <class T, class Policy>
 209   inline T hypergeometric_1F2_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol)
 210   {
 211     detail::hypergeometric_pFq_cf_term<T, 1u, 2u> f(a, b, c, z);
 212     T result = detail::compute_cf_pFq(f, pol);
 213     result = (((a * z) / (b * c)) / result) + 1;
 214     return result;
 215   }
 216
 217   template <class T, class Policy>
 218   inline T hypergeometric_2F1_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol)
 219   {
 220     detail::hypergeometric_pFq_cf_term<T, 2u, 1u> f(a, b, c, z);
 221     T result = detail::compute_cf_pFq(f, pol);
 222     result = ((((a * b) * z) / c) / result) + 1;
 223     return result;
 224   }
 225
 226   } } } // namespaces
 227
 228 #endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP