117 lines · plain
1///////////////////////////////////////////////////////////////////////////////2// Copyright 2014 Anton Bikineev3// Copyright 2014 Christopher Kormanyos4// Copyright 2014 John Maddock5// Copyright 2014 Paul Bristow6// Distributed under the Boost7// Software License, Version 1.0. (See accompanying file8// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)9 10#ifndef BOOST_MATH_HYPERGEOMETRIC_0F1_HPP11#define BOOST_MATH_HYPERGEOMETRIC_0F1_HPP12 13#include <boost/math/policies/policy.hpp>14#include <boost/math/policies/error_handling.hpp>15#include <boost/math/special_functions/detail/hypergeometric_series.hpp>16#include <boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp>17 18namespace boost { namespace math { namespace detail {19 20 21 template <class T>22 struct hypergeometric_0F1_cf23 {24 //25 // We start this continued fraction at b on index -126 // and treat the -1 and 0 cases as special cases.27 // We do this to avoid adding the continued fraction result28 // to 1 so that we can accurately evaluate for small results29 // as well as large ones. See http://functions.wolfram.com/07.17.10.0002.0130 //31 T b, z;32 int k;33 hypergeometric_0F1_cf(T b_, T z_) : b(b_), z(z_), k(-2) {}34 typedef std::pair<T, T> result_type;35 36 result_type operator()()37 {38 ++k;39 if (k <= 0)40 return std::make_pair(z / b, 1);41 return std::make_pair(-z / ((k + 1) * (b + k)), 1 + z / ((k + 1) * (b + k)));42 }43 };44 45 template <class T, class Policy>46 T hypergeometric_0F1_cf_imp(T b, T z, const Policy& pol, const char* function)47 {48 hypergeometric_0F1_cf<T> evaluator(b, z);49 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();50 T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter);51 policies::check_series_iterations<T>(function, max_iter, pol);52 return cf;53 }54 55 56 template <class T, class Policy>57 inline T hypergeometric_0F1_imp(const T& b, const T& z, const Policy& pol)58 {59 const char* function = "boost::math::hypergeometric_0f1<%1%,%1%>(%1%, %1%)";60 BOOST_MATH_STD_USING61 62 // some special cases63 if (z == 0)64 return T(1);65 66 if ((b <= 0) && (b == floor(b)))67 return policies::raise_pole_error<T>(function, "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol);68 69 if (z < -5 && b > -5)70 {71 // Series is alternating and divergent, need to do something else here,72 // Bessel function relation is much more accurate, unless |b| is similarly73 // large to |z|, otherwise the CF formula suffers from cancellation when74 // the result would be very small.75 if (fabs(z / b) > 4)76 return hypergeometric_0F1_bessel(b, z, pol);77 return hypergeometric_0F1_cf_imp(b, z, pol, function);78 }79 // evaluation through Taylor series looks80 // more precisious than Bessel relation:81 // detail::hypergeometric_0f1_bessel(b, z, pol);82 return detail::hypergeometric_0F1_generic_series(b, z, pol);83 }84 85} // namespace detail86 87template <class T1, class T2, class Policy>88inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z, const Policy& /* pol */)89{90 BOOST_FPU_EXCEPTION_GUARD91 typedef typename tools::promote_args<T1, T2>::type result_type;92 typedef typename policies::evaluation<result_type, Policy>::type value_type;93 typedef typename policies::normalise<94 Policy,95 policies::promote_float<false>,96 policies::promote_double<false>,97 policies::discrete_quantile<>,98 policies::assert_undefined<> >::type forwarding_policy;99 return policies::checked_narrowing_cast<result_type, Policy>(100 detail::hypergeometric_0F1_imp<value_type>(101 static_cast<value_type>(b),102 static_cast<value_type>(z),103 forwarding_policy()),104 "boost::math::hypergeometric_0F1<%1%>(%1%,%1%)");105}106 107template <class T1, class T2>108inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z)109{110 return hypergeometric_0F1(b, z, policies::policy<>());111}112 113 114} } // namespace boost::math115 116#endif // BOOST_MATH_HYPERGEOMETRIC_HPP117