164 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_2F0_HPP11#define BOOST_MATH_HYPERGEOMETRIC_2F0_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/laguerre.hpp>17#include <boost/math/special_functions/hermite.hpp>18#include <boost/math/tools/fraction.hpp>19 20namespace boost { namespace math { namespace detail {21 22 template <class T>23 struct hypergeometric_2F0_cf24 {25 //26 // We start this continued fraction at b on index -127 // and treat the -1 and 0 cases as special cases.28 // We do this to avoid adding the continued fraction result29 // to 1 so that we can accurately evaluate for small results30 // as well as large ones. See http://functions.wolfram.com/07.31.10.0002.0131 //32 T a1, a2, z;33 int k;34 hypergeometric_2F0_cf(T a1_, T a2_, T z_) : a1(a1_), a2(a2_), z(z_), k(-2) {}35 typedef std::pair<T, T> result_type;36 37 result_type operator()()38 {39 ++k;40 if (k <= 0)41 return std::make_pair(z * a1 * a2, 1);42 return std::make_pair(-z * (a1 + k) * (a2 + k) / (k + 1), 1 + z * (a1 + k) * (a2 + k) / (k + 1));43 }44 };45 46 template <class T, class Policy>47 T hypergeometric_2F0_cf_imp(T a1, T a2, T z, const Policy& pol, const char* function)48 {49 using namespace boost::math;50 hypergeometric_2F0_cf<T> evaluator(a1, a2, z);51 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();52 T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter);53 policies::check_series_iterations<T>(function, max_iter, pol);54 return cf;55 }56 57 58 template <class T, class Policy>59 inline T hypergeometric_2F0_imp(T a1, T a2, const T& z, const Policy& pol, bool asymptotic = false)60 {61 //62 // The terms in this series go to infinity unless one of a1 and a2 is a negative integer.63 //64 using std::swap;65 BOOST_MATH_STD_USING66 67 static const char* const function = "boost::math::hypergeometric_2F0<%1%,%1%,%1%>(%1%,%1%,%1%)";68 69 if (z == 0)70 return 1;71 72 bool is_a1_integer = (a1 == floor(a1));73 bool is_a2_integer = (a2 == floor(a2));74 75 if (!asymptotic && !is_a1_integer && !is_a2_integer)76 return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);77 if (!is_a1_integer || (a1 > 0))78 {79 swap(a1, a2);80 swap(is_a1_integer, is_a2_integer);81 }82 //83 // At this point a1 must be a negative integer:84 //85 if(!asymptotic && (!is_a1_integer || (a1 > 0)))86 return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);87 //88 // Special cases first:89 //90 if (a1 == 0)91 return 1;92 if ((a1 == a2 - 0.5f) && (z < 0))93 {94 // http://functions.wolfram.com/07.31.03.0083.0195 int n = static_cast<int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-2 * a1)));96 T smz = sqrt(-z);97 return static_cast<T>(pow(2 / smz, T(-n)) * boost::math::hermite(n, 1 / smz, pol)); // Warning suppression: integer power returns at least a double98 }99 100 if (is_a1_integer && is_a2_integer)101 {102 if ((a1 < 1) && (a2 <= a1))103 {104 const unsigned int n = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a1)));105 const unsigned int m = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a2 - n)));106 107 return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) *108 boost::math::laguerre(n, m, -(1 / z), pol);109 }110 else if ((a2 < 1) && (a1 <= a2))111 {112 // function is symmetric for a1 and a2113 const unsigned int n = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a2)));114 const unsigned int m = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a1 - n)));115 116 return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) *117 boost::math::laguerre(n, m, -(1 / z), pol);118 }119 }120 121 if ((a1 * a2 * z < 0) && (a2 < -5) && (fabs(a1 * a2 * z) > 0.5))122 {123 // Series is alternating and maybe divergent at least for the first few terms124 // (until a2 goes positive), try the continued fraction:125 return hypergeometric_2F0_cf_imp(a1, a2, z, pol, function);126 }127 128 return detail::hypergeometric_2F0_generic_series(a1, a2, z, pol);129 }130 131} // namespace detail132 133template <class T1, class T2, class T3, class Policy>134inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z, const Policy& /* pol */)135{136 BOOST_FPU_EXCEPTION_GUARD137 typedef typename tools::promote_args<T1, T2, T3>::type result_type;138 typedef typename policies::evaluation<result_type, Policy>::type value_type;139 typedef typename policies::normalise<140 Policy,141 policies::promote_float<false>,142 policies::promote_double<false>,143 policies::discrete_quantile<>,144 policies::assert_undefined<> >::type forwarding_policy;145 return policies::checked_narrowing_cast<result_type, Policy>(146 detail::hypergeometric_2F0_imp<value_type>(147 static_cast<value_type>(a1),148 static_cast<value_type>(a2),149 static_cast<value_type>(z),150 forwarding_policy()),151 "boost::math::hypergeometric_2F0<%1%>(%1%,%1%,%1%)");152}153 154template <class T1, class T2, class T3>155inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z)156{157 return hypergeometric_2F0(a1, a2, z, policies::policy<>());158}159 160 161 } } // namespace boost::math162 163#endif // BOOST_MATH_HYPERGEOMETRIC_HPP164