boost/math/special_functions/hypergeometric_2F0.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2014 Anton Bikineev
   3 //  Copyright 2014 Christopher Kormanyos
   4 //  Copyright 2014 John Maddock
   5 //  Copyright 2014 Paul Bristow
   6 //  Distributed under the Boost
   7 //  Software License, Version 1.0. (See accompanying file
   8 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   9
  10 #ifndef BOOST_MATH_HYPERGEOMETRIC_2F0_HPP
  11 #define BOOST_MATH_HYPERGEOMETRIC_2F0_HPP
  12
  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
  20 namespace boost { namespace math { namespace detail {
  21
  22    template <class T>
  23    struct hypergeometric_2F0_cf
  24    {
  25       //
  26       // We start this continued fraction at b on index -1
  27       // and treat the -1 and 0 cases as special cases.
  28       // We do this to avoid adding the continued fraction result
  29       // to 1 so that we can accurately evaluate for small results
  30       // as well as large ones.  See  http://functions.wolfram.com/07.31.10.0002.01
  31       //
  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       boost::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_USING
  66
  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, 0, 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, 0, 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.01
  95          int n = static_cast<int>(static_cast<boost::uintmax_t>(boost::math::lltrunc(-2 * a1)));
  96          T smz = sqrt(-z);
  97          return pow(2 / smz, -n) * boost::math::hermite(n, 1 / smz);
  98       }
  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<boost::uintmax_t>(boost::math::lltrunc(-a1)));
 105             const unsigned int m = static_cast<unsigned int>(static_cast<boost::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 a2
 113             const unsigned int n = static_cast<unsigned int>(static_cast<boost::uintmax_t>(boost::math::lltrunc(-a2)));
 114             const unsigned int m = static_cast<unsigned int>(static_cast<boost::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 terms
 124          // (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 detail
 132
 133 template <class T1, class T2, class T3, class Policy>
 134 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z, const Policy& /* pol */)
 135 {
 136    BOOST_FPU_EXCEPTION_GUARD
 137       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
 154 template <class T1, class T2, class T3>
 155 inline 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::math
 162
 163 #endif // BOOST_MATH_HYPERGEOMETRIC_HPP