math/special_functions/hypergeometric_2F0.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2014 Anton Bikineev
0003 //  Copyright 2014 Christopher Kormanyos
0004 //  Copyright 2014 John Maddock
0005 //  Copyright 2014 Paul Bristow
0006 //  Distributed under the Boost
0007 //  Software License, Version 1.0. (See accompanying file
0008 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0009
0010 #ifndef BOOST_MATH_HYPERGEOMETRIC_2F0_HPP
0011 #define BOOST_MATH_HYPERGEOMETRIC_2F0_HPP
0012
0013 #include <boost/math/policies/policy.hpp>
0014 #include <boost/math/policies/error_handling.hpp>
0015 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
0016 #include <boost/math/special_functions/laguerre.hpp>
0017 #include <boost/math/special_functions/hermite.hpp>
0018 #include <boost/math/tools/fraction.hpp>
0019
0020 namespace boost { namespace math { namespace detail {
0021
0022    template <class T>
0023    struct hypergeometric_2F0_cf
0024    {
0025       //
0026       // We start this continued fraction at b on index -1
0027       // and treat the -1 and 0 cases as special cases.
0028       // We do this to avoid adding the continued fraction result
0029       // to 1 so that we can accurately evaluate for small results
0030       // as well as large ones.  See  http://functions.wolfram.com/07.31.10.0002.01
0031       //
0032       T a1, a2, z;
0033       int k;
0034       hypergeometric_2F0_cf(T a1_, T a2_, T z_) : a1(a1_), a2(a2_), z(z_), k(-2) {}
0035       typedef std::pair<T, T> result_type;
0036
0037       result_type operator()()
0038       {
0039          ++k;
0040          if (k <= 0)
0041             return std::make_pair(z * a1 * a2, 1);
0042          return std::make_pair(-z * (a1 + k) * (a2 + k) / (k + 1), 1 + z * (a1 + k) * (a2 + k) / (k + 1));
0043       }
0044    };
0045
0046    template <class T, class Policy>
0047    T hypergeometric_2F0_cf_imp(T a1, T a2, T z, const Policy& pol, const char* function)
0048    {
0049       using namespace boost::math;
0050       hypergeometric_2F0_cf<T> evaluator(a1, a2, z);
0051       std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0052       T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter);
0053       policies::check_series_iterations<T>(function, max_iter, pol);
0054       return cf;
0055    }
0056
0057
0058    template <class T, class Policy>
0059    inline T hypergeometric_2F0_imp(T a1, T a2, const T& z, const Policy& pol, bool asymptotic = false)
0060    {
0061       //
0062       // The terms in this series go to infinity unless one of a1 and a2 is a negative integer.
0063       //
0064       using std::swap;
0065       BOOST_MATH_STD_USING
0066
0067       static const char* const function = "boost::math::hypergeometric_2F0<%1%,%1%,%1%>(%1%,%1%,%1%)";
0068
0069       if (z == 0)
0070          return 1;
0071
0072       bool is_a1_integer = (a1 == floor(a1));
0073       bool is_a2_integer = (a2 == floor(a2));
0074
0075       if (!asymptotic && !is_a1_integer && !is_a2_integer)
0076          return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
0077       if (!is_a1_integer || (a1 > 0))
0078       {
0079          swap(a1, a2);
0080          swap(is_a1_integer, is_a2_integer);
0081       }
0082       //
0083       // At this point a1 must be a negative integer:
0084       //
0085       if(!asymptotic && (!is_a1_integer || (a1 > 0)))
0086          return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
0087       //
0088       // Special cases first:
0089       //
0090       if (a1 == 0)
0091          return 1;
0092       if ((a1 == a2 - 0.5f) && (z < 0))
0093       {
0094          // http://functions.wolfram.com/07.31.03.0083.01
0095          int n = static_cast<int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-2 * a1)));
0096          T smz = sqrt(-z);
0097          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
0098       }
0099
0100       if (is_a1_integer && is_a2_integer)
0101       {
0102          if ((a1 < 1) && (a2 <= a1))
0103          {
0104             const unsigned int n = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a1)));
0105             const unsigned int m = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a2 - n)));
0106
0107             return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) *
0108                boost::math::laguerre(n, m, -(1 / z), pol);
0109          }
0110          else if ((a2 < 1) && (a1 <= a2))
0111          {
0112             // function is symmetric for a1 and a2
0113             const unsigned int n = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a2)));
0114             const unsigned int m = static_cast<unsigned int>(static_cast<std::uintmax_t>(boost::math::lltrunc(-a1 - n)));
0115
0116             return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) *
0117                boost::math::laguerre(n, m, -(1 / z), pol);
0118          }
0119       }
0120
0121       if ((a1 * a2 * z < 0) && (a2 < -5) && (fabs(a1 * a2 * z) > 0.5))
0122       {
0123          // Series is alternating and maybe divergent at least for the first few terms
0124          // (until a2 goes positive), try the continued fraction:
0125          return hypergeometric_2F0_cf_imp(a1, a2, z, pol, function);
0126       }
0127
0128       return detail::hypergeometric_2F0_generic_series(a1, a2, z, pol);
0129    }
0130
0131 } // namespace detail
0132
0133 template <class T1, class T2, class T3, class Policy>
0134 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z, const Policy& /* pol */)
0135 {
0136    BOOST_FPU_EXCEPTION_GUARD
0137       typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0138    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0139    typedef typename policies::normalise<
0140       Policy,
0141       policies::promote_float<false>,
0142       policies::promote_double<false>,
0143       policies::discrete_quantile<>,
0144       policies::assert_undefined<> >::type forwarding_policy;
0145    return policies::checked_narrowing_cast<result_type, Policy>(
0146       detail::hypergeometric_2F0_imp<value_type>(
0147          static_cast<value_type>(a1),
0148          static_cast<value_type>(a2),
0149          static_cast<value_type>(z),
0150          forwarding_policy()),
0151       "boost::math::hypergeometric_2F0<%1%>(%1%,%1%,%1%)");
0152 }
0153
0154 template <class T1, class T2, class T3>
0155 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z)
0156 {
0157    return hypergeometric_2F0(a1, a2, z, policies::policy<>());
0158 }
0159
0160
0161   } } // namespace boost::math
0162
0163 #endif // BOOST_MATH_HYPERGEOMETRIC_HPP