boost/math/special_functions/detail/hypergeometric_0F1_bessel.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_0F1_BESSEL_HPP
  11 #define BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP
  12
  13 #include <boost/math/special_functions/bessel.hpp>
  14 #include <boost/math/special_functions/gamma.hpp>
  15
  16   namespace boost { namespace math { namespace detail {
  17
  18   template <class T, class Policy>
  19   inline T hypergeometric_0F1_bessel(const T& b, const T& z, const Policy& pol)
  20   {
  21     BOOST_MATH_STD_USING
  22
  23     const bool is_z_nonpositive = z <= 0;
  24
  25     const T sqrt_z = is_z_nonpositive ? T(sqrt(-z)) : T(sqrt(z));
  26     const T bessel_mult = is_z_nonpositive ?
  27       boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol) :
  28       boost::math::cyl_bessel_i(b - 1, 2 * sqrt_z, pol) ;
  29
  30     if (b > boost::math::max_factorial<T>::value)
  31     {
  32        const T lsqrt_z = log(sqrt_z);
  33        const T lsqrt_z_pow_b = (b - 1) * lsqrt_z;
  34        T lg = (boost::math::lgamma(b, pol) - lsqrt_z_pow_b);
  35        lg = exp(lg);
  36        return lg * bessel_mult;
  37     }
  38     else
  39     {
  40        const T sqrt_z_pow_b = pow(sqrt_z, b - 1);
  41        return (boost::math::tgamma(b, pol) / sqrt_z_pow_b) * bessel_mult;
  42     }
  43   }
  44
  45   } } } // namespaces
  46
  47 #endif // BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP