boost/math/special_functions/detail/hypergeometric_1F1_small_a_negative_b_by_ratio.hpp

   1
   2 ///////////////////////////////////////////////////////////////////////////////
   3 //  Copyright 2018 John Maddock
   4 //  Distributed under the Boost
   5 //  Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7 //
   8 #ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_SMALL_A_NEG_B_HPP
   9 #define BOOST_MATH_HYPERGEOMETRIC_1F1_SMALL_A_NEG_B_HPP
  10
  11 #include <algorithm>
  12 #include <boost/math/tools/recurrence.hpp>
  13
  14   namespace boost { namespace math { namespace detail {
  15
  16      // forward declaration for initial values
  17      template <class T, class Policy>
  18      inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol);
  19
  20      template <class T, class Policy>
  21      inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling);
  22
  23       template <class T>
  24       T max_b_for_1F1_small_a_negative_b_by_ratio(const T& z)
  25       {
  26          if (z < -998)
  27             return (z * 2) / 3;
  28          float max_b[][2] =
  29          {
  30             { 0.0f, -47.3046f }, {-6.7275f, -52.0351f }, { -8.9543f, -57.2386f }, {-11.9182f, -62.9625f }, {-14.421f, -69.2587f }, {-19.1943f, -76.1846f }, {-23.2252f, -83.803f }, {-28.1024f, -92.1833f }, {-34.0039f, -101.402f }, {-37.4043f, -111.542f }, {-45.2593f, -122.696f }, {-54.7637f, -134.966f }, {-60.2401f, -148.462f }, {-72.8905f, -163.308f }, {-88.1975f, -179.639f }, {-88.1975f, -197.603f }, {-106.719f, -217.363f }, {-129.13f, -239.1f }, {-142.043f, -263.01f }, {-156.247f, -289.311f }, {-189.059f, -318.242f }, {-207.965f, -350.066f }, {-228.762f, -385.073f }, {-276.801f, -423.58f }, {-304.482f, -465.938f }, {-334.93f, -512.532f }, {-368.423f, -563.785f }, {-405.265f, -620.163f }, {-445.792f, -682.18f }, {-539.408f, -750.398f }, {-593.349f, -825.437f }, {-652.683f, -907.981f }, {-717.952f, -998.779f }
  31          };
  32          auto p = std::lower_bound(max_b, max_b + sizeof(max_b) / sizeof(max_b[0]), z, [](const float (&a)[2], const T& z) { return a[1] > z; });
  33          T b = p - max_b ? (*--p)[0] : 0;
  34          //
  35          // We need approximately an extra 10 recurrences per 50 binary digits precision above that of double:
  36          //
  37          b += (std::max)(0, boost::math::tools::digits<T>() - 53) / 5;
  38          return b;
  39       }
  40
  41       template <class T, class Policy>
  42       T hypergeometric_1F1_small_a_negative_b_by_ratio(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
  43       {
  44          BOOST_MATH_STD_USING
  45          //
  46          // We grab the ratio for M[a, b, z] / M[a, b+1, z] and use it to seed 2 initial values,
  47          // then recurse until b > 0, compute a reference value and normalize (Millers method).
  48          //
  49          int iterations = itrunc(-b, pol);
  50          boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
  51          T ratio = boost::math::tools::function_ratio_from_forwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T>(a, b, z), boost::math::tools::epsilon<T>(), max_iter);
  52          boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_small_a_negative_b_by_ratio<%1%>(%1%,%1%,%1%)", max_iter, pol);
  53          T first = 1;
  54          T second = 1 / ratio;
  55          int scaling1 = 0;
  56          BOOST_ASSERT(b + iterations != a);
  57          second = boost::math::tools::apply_recurrence_relation_forward(boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T>(a, b + 1, z), iterations, first, second, &scaling1);
  58          int scaling2 = 0;
  59          first = hypergeometric_1F1_imp(a, T(b + iterations + 1), z, pol, scaling2);
  60          //
  61          // Result is now first/second * e^(scaling2 - scaling1)
  62          //
  63          log_scaling += scaling2 - scaling1;
  64          return first / second;
  65       }
  66
  67
  68   } } } // namespaces
  69
  70 #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_SMALL_A_NEG_B_HPP