special_functions/detail/hypergeometric_cf.hpp

0001
0002 ///////////////////////////////////////////////////////////////////////////////
0003 //  Copyright 2014 Anton Bikineev
0004 //  Copyright 2014 Christopher Kormanyos
0005 //  Copyright 2014 John Maddock
0006 //  Copyright 2014 Paul Bristow
0007 //  Distributed under the Boost
0008 //  Software License, Version 1.0. (See accompanying file
0009 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0010 //
0011 #ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP
0012 #define BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP
0013
0014   namespace boost { namespace math { namespace detail {
0015
0016   // primary template for term of continued fraction
0017   template <class T, unsigned p, unsigned q>
0018   struct hypergeometric_pFq_cf_term;
0019
0020   // partial specialization for 0F1
0021   template <class T>
0022   struct hypergeometric_pFq_cf_term<T, 0u, 1u>
0023   {
0024     typedef std::pair<T,T> result_type;
0025
0026     hypergeometric_pFq_cf_term(const T& b, const T& z):
0027       n(1), b(b), z(z),
0028       term(std::make_pair(T(0), T(1)))
0029     {
0030     }
0031
0032     result_type operator()()
0033     {
0034       const result_type result = term;
0035       ++b; ++n;
0036       numer = -(z / (b * n));
0037       term = std::make_pair(numer, 1 - numer);
0038       return result;
0039     }
0040
0041   private:
0042     unsigned n;
0043     T b;
0044     const T z;
0045     T numer;
0046     result_type term;
0047   };
0048
0049   // partial specialization for 1F0
0050   template <class T>
0051   struct hypergeometric_pFq_cf_term<T, 1u, 0u>
0052   {
0053     typedef std::pair<T,T> result_type;
0054
0055     hypergeometric_pFq_cf_term(const T& a, const T& z):
0056       n(1), a(a), z(z),
0057       term(std::make_pair(T(0), T(1)))
0058     {
0059     }
0060
0061     result_type operator()()
0062     {
0063       const result_type result = term;
0064       ++a; ++n;
0065       numer = -((a * z) / n);
0066       term = std::make_pair(numer, 1 - numer);
0067       return result;
0068     }
0069
0070   private:
0071     unsigned n;
0072     T a;
0073     const T z;
0074     T numer;
0075     result_type term;
0076   };
0077
0078   // partial specialization for 1F1
0079   template <class T>
0080   struct hypergeometric_pFq_cf_term<T, 1u, 1u>
0081   {
0082     typedef std::pair<T,T> result_type;
0083
0084     hypergeometric_pFq_cf_term(const T& a, const T& b, const T& z):
0085       n(1), a(a), b(b), z(z),
0086       term(std::make_pair(T(0), T(1)))
0087     {
0088     }
0089
0090     result_type operator()()
0091     {
0092       const result_type result = term;
0093       ++a; ++b; ++n;
0094       numer = -((a * z) / (b * n));
0095       term = std::make_pair(numer, 1 - numer);
0096       return result;
0097     }
0098
0099   private:
0100     unsigned n;
0101     T a, b;
0102     const T z;
0103     T numer;
0104     result_type term;
0105   };
0106
0107   // partial specialization for 1f2
0108   template <class T>
0109   struct hypergeometric_pFq_cf_term<T, 1u, 2u>
0110   {
0111     typedef std::pair<T,T> result_type;
0112
0113     hypergeometric_pFq_cf_term(const T& a, const T& b, const T& c, const T& z):
0114       n(1), a(a), b(b), c(c), z(z),
0115       term(std::make_pair(T(0), T(1)))
0116     {
0117     }
0118
0119     result_type operator()()
0120     {
0121       const result_type result = term;
0122       ++a; ++b; ++c; ++n;
0123       numer = -((a * z) / ((b * c) * n));
0124       term = std::make_pair(numer, 1 - numer);
0125       return result;
0126     }
0127
0128   private:
0129     unsigned n;
0130     T a, b, c;
0131     const T z;
0132     T numer;
0133     result_type term;
0134   };
0135
0136   // partial specialization for 2f1
0137   template <class T>
0138   struct hypergeometric_pFq_cf_term<T, 2u, 1u>
0139   {
0140     typedef std::pair<T,T> result_type;
0141
0142     hypergeometric_pFq_cf_term(const T& a, const T& b, const T& c, const T& z):
0143       n(1), a(a), b(b), c(c), z(z),
0144       term(std::make_pair(T(0), T(1)))
0145     {
0146     }
0147
0148     result_type operator()()
0149     {
0150       const result_type result = term;
0151       ++a; ++b; ++c; ++n;
0152       numer = -(((a * b) * z) / (c * n));
0153       term = std::make_pair(numer, 1 - numer);
0154       return result;
0155     }
0156
0157   private:
0158     unsigned n;
0159     T a, b, c;
0160     const T z;
0161     T numer;
0162     result_type term;
0163   };
0164
0165   template <class T, unsigned p, unsigned q, class Policy>
0166   inline T compute_cf_pFq(detail::hypergeometric_pFq_cf_term<T, p, q>& term, const Policy& pol)
0167   {
0168     BOOST_MATH_STD_USING
0169     std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0170     const T result = tools::continued_fraction_b(
0171       term,
0172       boost::math::policies::get_epsilon<T, Policy>(),
0173       max_iter);
0174     boost::math::policies::check_series_iterations<T>(
0175       "boost::math::hypergeometric_pFq_cf<%1%>(%1%,%1%,%1%)",
0176       max_iter,
0177       pol);
0178     return result;
0179   }
0180
0181   template <class T, class Policy>
0182   inline T hypergeometric_0F1_cf(const T& b, const T& z, const Policy& pol)
0183   {
0184     detail::hypergeometric_pFq_cf_term<T, 0u, 1u> f(b, z);
0185     T result = detail::compute_cf_pFq(f, pol);
0186     result = ((z / b) / result) + 1;
0187     return result;
0188   }
0189
0190   template <class T, class Policy>
0191   inline T hypergeometric_1F0_cf(const T& a, const T& z, const Policy& pol)
0192   {
0193     detail::hypergeometric_pFq_cf_term<T, 1u, 0u> f(a, z);
0194     T result = detail::compute_cf_pFq(f, pol);
0195     result = ((a * z) / result) + 1;
0196     return result;
0197   }
0198
0199   template <class T, class Policy>
0200   inline T hypergeometric_1F1_cf(const T& a, const T& b, const T& z, const Policy& pol)
0201   {
0202     detail::hypergeometric_pFq_cf_term<T, 1u, 1u> f(a, b, z);
0203     T result = detail::compute_cf_pFq(f, pol);
0204     result = (((a * z) / b) / result) + 1;
0205     return result;
0206   }
0207
0208   template <class T, class Policy>
0209   inline T hypergeometric_1F2_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol)
0210   {
0211     detail::hypergeometric_pFq_cf_term<T, 1u, 2u> f(a, b, c, z);
0212     T result = detail::compute_cf_pFq(f, pol);
0213     result = (((a * z) / (b * c)) / result) + 1;
0214     return result;
0215   }
0216
0217   template <class T, class Policy>
0218   inline T hypergeometric_2F1_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol)
0219   {
0220     detail::hypergeometric_pFq_cf_term<T, 2u, 1u> f(a, b, c, z);
0221     T result = detail::compute_cf_pFq(f, pol);
0222     result = ((((a * b) * z) / c) / result) + 1;
0223     return result;
0224   }
0225
0226   } } } // namespaces
0227
0228 #endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP