/usr/include/boost/math/special_functions
Edit: /usr/include/boost/math/special_functions/hypergeometric_2F0.hpp (6127B)
///////////////////////////////////////////////////////////////////////////////
// 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
#include
#include
#include
#include
#include
namespace boost { namespace math { namespace detail {
template
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 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
T hypergeometric_2F0_cf_imp(T a1, T a2, T z, const Policy& pol, const char* function)
{
using namespace boost::math;
hypergeometric_2F0_cf evaluator(a1, a2, z);
boost::uintmax_t max_iter = policies::get_max_series_iterations();
T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon(), max_iter);
policies::check_series_iterations(function, max_iter, pol);
return cf;
}
template
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(function, 0, 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(function, 0, 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(static_cast(boost::math::lltrunc(-2 * a1)));
T smz = sqrt(-z);
return pow(2 / smz, -n) * boost::math::hermite(n, 1 / smz);
}
if (is_a1_integer && is_a2_integer)
{
if ((a1 < 1) && (a2 <= a1))
{
const unsigned int n = static_cast(static_cast(boost::math::lltrunc(-a1)));
const unsigned int m = static_cast(static_cast(boost::math::lltrunc(-a2 - n)));
return (pow(z, T(n)) * boost::math::factorial(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(static_cast(boost::math::lltrunc(-a2)));
const unsigned int m = static_cast(static_cast(boost::math::lltrunc(-a1 - n)));
return (pow(z, T(n)) * boost::math::factorial(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
inline typename tools::promote_args::type hypergeometric_2F0(T1 a1, T2 a2, T3 z, const Policy& /* pol */)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args::type result_type;
typedef typename policies::evaluation::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float,
policies::promote_double,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast(
detail::hypergeometric_2F0_imp(
static_cast(a1),
static_cast(a2),
static_cast(z),
forwarding_policy()),
"boost::math::hypergeometric_2F0<%1%>(%1%,%1%,%1%)");
}
template
inline typename tools::promote_args::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