math/special_functions/hypergeometric_0F1.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_0F1_HPP
0011 #define BOOST_MATH_HYPERGEOMETRIC_0F1_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/detail/hypergeometric_0F1_bessel.hpp>
0017
0018 namespace boost { namespace math { namespace detail {
0019
0020
0021    template <class T>
0022    struct hypergeometric_0F1_cf
0023    {
0024       //
0025       // We start this continued fraction at b on index -1
0026       // and treat the -1 and 0 cases as special cases.
0027       // We do this to avoid adding the continued fraction result
0028       // to 1 so that we can accurately evaluate for small results
0029       // as well as large ones.  See http://functions.wolfram.com/07.17.10.0002.01
0030       //
0031       T b, z;
0032       int k;
0033       hypergeometric_0F1_cf(T b_, T z_) : b(b_), z(z_), k(-2) {}
0034       typedef std::pair<T, T> result_type;
0035
0036       result_type operator()()
0037       {
0038          ++k;
0039          if (k <= 0)
0040             return std::make_pair(z / b, 1);
0041          return std::make_pair(-z / ((k + 1) * (b + k)), 1 + z / ((k + 1) * (b + k)));
0042       }
0043    };
0044
0045    template <class T, class Policy>
0046    T hypergeometric_0F1_cf_imp(T b, T z, const Policy& pol, const char* function)
0047    {
0048       hypergeometric_0F1_cf<T> evaluator(b, z);
0049       std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0050       T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter);
0051       policies::check_series_iterations<T>(function, max_iter, pol);
0052       return cf;
0053    }
0054
0055
0056    template <class T, class Policy>
0057    inline T hypergeometric_0F1_imp(const T& b, const T& z, const Policy& pol)
0058    {
0059       const char* function = "boost::math::hypergeometric_0f1<%1%,%1%>(%1%, %1%)";
0060       BOOST_MATH_STD_USING
0061
0062          // some special cases
0063          if (z == 0)
0064             return T(1);
0065
0066       if ((b <= 0) && (b == floor(b)))
0067          return policies::raise_pole_error<T>(
0068             function,
0069             "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol);
0070
0071       if (z < -5 && b > -5)
0072       {
0073          // Series is alternating and divergent, need to do something else here,
0074          // Bessel function relation is much more accurate, unless |b| is similarly
0075          // large to |z|, otherwise the CF formula suffers from cancellation when
0076          // the result would be very small.
0077          if (fabs(z / b) > 4)
0078             return hypergeometric_0F1_bessel(b, z, pol);
0079          return hypergeometric_0F1_cf_imp(b, z, pol, function);
0080       }
0081       // evaluation through Taylor series looks
0082       // more precisious than Bessel relation:
0083       // detail::hypergeometric_0f1_bessel(b, z, pol);
0084       return detail::hypergeometric_0F1_generic_series(b, z, pol);
0085    }
0086
0087 } // namespace detail
0088
0089 template <class T1, class T2, class Policy>
0090 inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z, const Policy& /* pol */)
0091 {
0092    BOOST_FPU_EXCEPTION_GUARD
0093       typedef typename tools::promote_args<T1, T2>::type result_type;
0094    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0095    typedef typename policies::normalise<
0096       Policy,
0097       policies::promote_float<false>,
0098       policies::promote_double<false>,
0099       policies::discrete_quantile<>,
0100       policies::assert_undefined<> >::type forwarding_policy;
0101    return policies::checked_narrowing_cast<result_type, Policy>(
0102       detail::hypergeometric_0F1_imp<value_type>(
0103          static_cast<value_type>(b),
0104          static_cast<value_type>(z),
0105          forwarding_policy()),
0106       "boost::math::hypergeometric_0F1<%1%>(%1%,%1%)");
0107 }
0108
0109 template <class T1, class T2>
0110 inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z)
0111 {
0112    return hypergeometric_0F1(b, z, policies::policy<>());
0113 }
0114
0115
0116 } } // namespace boost::math
0117
0118 #endif // BOOST_MATH_HYPERGEOMETRIC_HPP