229 lines · plain
1 2///////////////////////////////////////////////////////////////////////////////3// Copyright 2014 Anton Bikineev4// Copyright 2014 Christopher Kormanyos5// Copyright 2014 John Maddock6// Copyright 2014 Paul Bristow7// Distributed under the Boost8// Software License, Version 1.0. (See accompanying file9// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)10//11#ifndef BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP12#define BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP13 14 namespace boost { namespace math { namespace detail {15 16 // primary template for term of continued fraction17 template <class T, unsigned p, unsigned q>18 struct hypergeometric_pFq_cf_term;19 20 // partial specialization for 0F121 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 1F050 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 1F179 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 1f2108 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 2f1137 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_USING169 std::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 } } } // namespaces227 228#endif // BOOST_MATH_DETAIL_HYPERGEOMETRIC_CF_HPP229