hypergeometric_1F1_cf.hpp 1.8 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // Copyright 2018 John Maddock
  3. // Distributed under the Boost
  4. // Software License, Version 1.0. (See accompanying file
  5. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. //
  7. #ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_CF_HPP
  8. #define BOOST_MATH_HYPERGEOMETRIC_1F1_CF_HPP
  9. #include <boost/math/tools/fraction.hpp>
  10. //
  11. // Evaluation of 1F1 by continued fraction
  12. // see http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/10/0002/
  13. //
  14. // This is not terribly useful, as like the series we're adding a something to 1,
  15. // so only really useful when we know that the result will be > 1.
  16. //
  17. namespace boost { namespace math { namespace detail {
  18. template <class T>
  19. struct hypergeometric_1F1_cf_func
  20. {
  21. typedef std::pair<T, T> result_type;
  22. hypergeometric_1F1_cf_func(T a_, T b_, T z_) : a(a_), b(b_), z(z_), k(0) {}
  23. std::pair<T, T> operator()()
  24. {
  25. ++k;
  26. return std::make_pair(-(((a + k) * z) / ((k + 1) * (b + k))), 1 + ((a + k) * z) / ((k + 1) * (b + k)));
  27. }
  28. T a, b, z;
  29. unsigned k;
  30. };
  31. template <class T, class Policy>
  32. T hypergeometric_1F1_cf(const T& a, const T& b, const T& z, const Policy& pol, const char* function)
  33. {
  34. hypergeometric_1F1_cf_func<T> func(a, b, z);
  35. std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
  36. T result = boost::math::tools::continued_fraction_a(func, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
  37. boost::math::policies::check_series_iterations<T>(function, max_iter, pol);
  38. return 1 + a * z / (b * (1 + result));
  39. }
  40. } } } // namespaces
  41. #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP