special_functions/detail/hypergeometric_0F1_bessel.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_BESSEL_HPP
0011 #define BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP
0012
0013 #include <boost/math/special_functions/bessel.hpp>
0014 #include <boost/math/special_functions/gamma.hpp>
0015
0016   namespace boost { namespace math { namespace detail {
0017
0018   template <class T, class Policy>
0019   inline T hypergeometric_0F1_bessel(const T& b, const T& z, const Policy& pol)
0020   {
0021     BOOST_MATH_STD_USING
0022
0023     const bool is_z_nonpositive = z <= 0;
0024
0025     const T sqrt_z = is_z_nonpositive ? T(sqrt(-z)) : T(sqrt(z));
0026     const T bessel_mult = is_z_nonpositive ?
0027       boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol) :
0028       boost::math::cyl_bessel_i(b - 1, 2 * sqrt_z, pol) ;
0029
0030     if (b > boost::math::max_factorial<T>::value)
0031     {
0032        const T lsqrt_z = log(sqrt_z);
0033        const T lsqrt_z_pow_b = (b - 1) * lsqrt_z;
0034        T lg = (boost::math::lgamma(b, pol) - lsqrt_z_pow_b);
0035        lg = exp(lg);
0036        return lg * bessel_mult;
0037     }
0038     else
0039     {
0040        const T sqrt_z_pow_b = pow(sqrt_z, b - 1);
0041        return (boost::math::tgamma(b, pol) / sqrt_z_pow_b) * bessel_mult;
0042     }
0043   }
0044
0045   } } } // namespaces
0046
0047 #endif // BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP