boost/math/special_functions/hypergeometric_2F0.hpp
/////////////////////////////////////////////////////////////////////////////// // Copyright 2014 Anton Bikineev // Copyright 2014 Christopher Kormanyos // Copyright 2014 John Maddock // Copyright 2014 Paul Bristow // Distributed under the Boost // Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_MATH_HYPERGEOMETRIC_2F0_HPP #define BOOST_MATH_HYPERGEOMETRIC_2F0_HPP #include <boost/math/policies/policy.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/special_functions/detail/hypergeometric_series.hpp> #include <boost/math/special_functions/laguerre.hpp> #include <boost/math/special_functions/hermite.hpp> #include <boost/math/tools/fraction.hpp> namespace boost { namespace math { namespace detail { template <class T> struct hypergeometric_2F0_cf { // // We start this continued fraction at b on index -1 // and treat the -1 and 0 cases as special cases. // We do this to avoid adding the continued fraction result // to 1 so that we can accurately evaluate for small results // as well as large ones. See http://functions.wolfram.com/07.31.10.0002.01 // T a1, a2, z; int k; hypergeometric_2F0_cf(T a1_, T a2_, T z_) : a1(a1_), a2(a2_), z(z_), k(-2) {} typedef std::pair<T, T> result_type; result_type operator()() { ++k; if (k <= 0) return std::make_pair(z * a1 * a2, 1); return std::make_pair(-z * (a1 + k) * (a2 + k) / (k + 1), 1 + z * (a1 + k) * (a2 + k) / (k + 1)); } }; template <class T, class Policy> T hypergeometric_2F0_cf_imp(T a1, T a2, T z, const Policy& pol, const char* function) { using namespace boost::math; hypergeometric_2F0_cf<T> evaluator(a1, a2, z); std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter); policies::check_series_iterations<T>(function, max_iter, pol); return cf; } template <class T, class Policy> inline T hypergeometric_2F0_imp(T a1, T a2, const T& z, const Policy& pol, bool asymptotic = false) { // // The terms in this series go to infinity unless one of a1 and a2 is a negative integer. // using std::swap; BOOST_MATH_STD_USING static const char* const function = "boost::math::hypergeometric_2F0<%1%,%1%,%1%>(%1%,%1%,%1%)"; if (z == 0) return 1; bool is_a1_integer = (a1 == floor(a1)); bool is_a2_integer = (a2 == floor(a2)); if (!asymptotic && !is_a1_integer && !is_a2_integer) return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol); if (!is_a1_integer || (a1 > 0)) { swap(a1, a2); swap(is_a1_integer, is_a2_integer); } // // At this point a1 must be a negative integer: // if(!asymptotic && (!is_a1_integer || (a1 > 0))) return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol); // // Special cases first: // if (a1 == 0) return 1; if ((a1 == a2 - 0.5f) && (z < 0)) { // http://functions.wolfram.com/07.31.03.0083.01 int n = static_cast<int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-2 * a1))); T smz = sqrt(-z); return static_cast<T>(pow(2 / smz, T(-n)) * boost::math::hermite(n, 1 / smz, pol)); // Warning suppression: integer power returns at least a double } if (is_a1_integer && is_a2_integer) { if ((a1 < 1) && (a2 <= a1)) { const unsigned int n = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a1))); const unsigned int m = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a2 - n))); return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) * boost::math::laguerre(n, m, -(1 / z), pol); } else if ((a2 < 1) && (a1 <= a2)) { // function is symmetric for a1 and a2 const unsigned int n = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a2))); const unsigned int m = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a1 - n))); return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) * boost::math::laguerre(n, m, -(1 / z), pol); } } if ((a1 * a2 * z < 0) && (a2 < -5) && (fabs(a1 * a2 * z) > 0.5)) { // Series is alternating and maybe divergent at least for the first few terms // (until a2 goes positive), try the continued fraction: return hypergeometric_2F0_cf_imp(a1, a2, z, pol, function); } return detail::hypergeometric_2F0_generic_series(a1, a2, z, pol); } } // namespace detail template <class T1, class T2, class T3, class Policy> inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z, const Policy& /* pol */) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args<T1, T2, T3>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; typedef typename policies::normalise< Policy, policies::promote_float<false>, policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; return policies::checked_narrowing_cast<result_type, Policy>( detail::hypergeometric_2F0_imp<value_type>( static_cast<value_type>(a1), static_cast<value_type>(a2), static_cast<value_type>(z), forwarding_policy()), "boost::math::hypergeometric_2F0<%1%>(%1%,%1%,%1%)"); } template <class T1, class T2, class T3> inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z) { return hypergeometric_2F0(a1, a2, z, policies::policy<>()); } } } // namespace boost::math #endif // BOOST_MATH_HYPERGEOMETRIC_HPP