hypergeometric_0F1.hpp 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118
  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. #ifndef BOOST_MATH_HYPERGEOMETRIC_0F1_HPP
  10. #define BOOST_MATH_HYPERGEOMETRIC_0F1_HPP
  11. #include <boost/math/policies/policy.hpp>
  12. #include <boost/math/policies/error_handling.hpp>
  13. #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
  14. #include <boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp>
  15. namespace boost { namespace math { namespace detail {
  16. template <class T>
  17. struct hypergeometric_0F1_cf
  18. {
  19. //
  20. // We start this continued fraction at b on index -1
  21. // and treat the -1 and 0 cases as special cases.
  22. // We do this to avoid adding the continued fraction result
  23. // to 1 so that we can accurately evaluate for small results
  24. // as well as large ones. See http://functions.wolfram.com/07.17.10.0002.01
  25. //
  26. T b, z;
  27. int k;
  28. hypergeometric_0F1_cf(T b_, T z_) : b(b_), z(z_), k(-2) {}
  29. typedef std::pair<T, T> result_type;
  30. result_type operator()()
  31. {
  32. ++k;
  33. if (k <= 0)
  34. return std::make_pair(z / b, 1);
  35. return std::make_pair(-z / ((k + 1) * (b + k)), 1 + z / ((k + 1) * (b + k)));
  36. }
  37. };
  38. template <class T, class Policy>
  39. T hypergeometric_0F1_cf_imp(T b, T z, const Policy& pol, const char* function)
  40. {
  41. hypergeometric_0F1_cf<T> evaluator(b, z);
  42. boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
  43. T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter);
  44. policies::check_series_iterations<T>(function, max_iter, pol);
  45. return cf;
  46. }
  47. template <class T, class Policy>
  48. inline T hypergeometric_0F1_imp(const T& b, const T& z, const Policy& pol)
  49. {
  50. const char* function = "boost::math::hypergeometric_0f1<%1%,%1%>(%1%, %1%)";
  51. BOOST_MATH_STD_USING
  52. // some special cases
  53. if (z == 0)
  54. return T(1);
  55. if ((b <= 0) && (b == floor(b)))
  56. return policies::raise_pole_error<T>(
  57. function,
  58. "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol);
  59. if (z < -5 && b > -5)
  60. {
  61. // Series is alternating and divergent, need to do something else here,
  62. // Bessel function relation is much more accurate, unless |b| is similarly
  63. // large to |z|, otherwise the CF formula suffers from cancellation when
  64. // the result would be very small.
  65. if (fabs(z / b) > 4)
  66. return hypergeometric_0F1_bessel(b, z, pol);
  67. return hypergeometric_0F1_cf_imp(b, z, pol, function);
  68. }
  69. // evaluation through Taylor series looks
  70. // more precisious than Bessel relation:
  71. // detail::hypergeometric_0f1_bessel(b, z, pol);
  72. return detail::hypergeometric_0F1_generic_series(b, z, pol);
  73. }
  74. } // namespace detail
  75. template <class T1, class T2, class Policy>
  76. inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z, const Policy& /* pol */)
  77. {
  78. BOOST_FPU_EXCEPTION_GUARD
  79. typedef typename tools::promote_args<T1, T2>::type result_type;
  80. typedef typename policies::evaluation<result_type, Policy>::type value_type;
  81. typedef typename policies::normalise<
  82. Policy,
  83. policies::promote_float<false>,
  84. policies::promote_double<false>,
  85. policies::discrete_quantile<>,
  86. policies::assert_undefined<> >::type forwarding_policy;
  87. return policies::checked_narrowing_cast<result_type, Policy>(
  88. detail::hypergeometric_0F1_imp<value_type>(
  89. static_cast<value_type>(b),
  90. static_cast<value_type>(z),
  91. forwarding_policy()),
  92. "boost::math::hypergeometric_0F1<%1%>(%1%,%1%)");
  93. }
  94. template <class T1, class T2>
  95. inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z)
  96. {
  97. return hypergeometric_0F1(b, z, policies::policy<>());
  98. }
  99. } } // namespace boost::math
  100. #endif // BOOST_MATH_HYPERGEOMETRIC_HPP