From a187b714e95abef047ccd628853fa64786b7987c Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Thu, 31 Jan 2019 19:06:17 +0000
Subject: [PATCH] 1F1: Apply backwards recurrence relations for GammaP in large
 a,b,z approximation. [CI SKIP]

---
 .../detail/hypergeometric_1F1_large_abz.hpp   | 24 +++++++++++++++++--
 .../detail/hypergeometric_1F1_recurrence.hpp  |  1 -
 test/test_1F1.hpp                             |  6 ++++-
 3 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/include/boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp b/include/boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp
index 9b54f2a51..e8b4d7c1d 100644
--- a/include/boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp
+++ b/include/boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp
@@ -19,26 +19,46 @@
      template <class T, class Policy>
      struct hypergeometric_1F1_igamma_series
      {
+        enum{ cache_size = 64 };
+
         typedef T result_type;
         hypergeometric_1F1_igamma_series(const T& alpha, const T& delta, const T& x, const Policy& pol)
-           : delta_poch(-delta), alpha_poch(alpha), x(x), k(0), pol(pol)
+           : delta_poch(-delta), alpha_poch(alpha), x(x), k(0), cache_offset(0), pol(pol)
         {
            BOOST_MATH_STD_USING
            T log_term = log(x) * -alpha;
            log_scaling = itrunc(log_term - 3 - boost::math::tools::log_min_value<T>() / 50);
            term = exp(log_term - log_scaling);
+           refill_cache();
         }
         T operator()()
         {
-           T result = term * boost::math::gamma_p(alpha_poch, x, pol);
+           if (k - cache_offset >= cache_size)
+           {
+              cache_offset += cache_size;
+              refill_cache();
+           }
+           T result = term * gamma_cache[k - cache_offset];
            term *= delta_poch * alpha_poch / (++k * x);
            delta_poch += 1;
            alpha_poch += 1;
            return result;
         }
+        void refill_cache()
+        {
+           typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
+
+           gamma_cache[cache_size - 1] = boost::math::gamma_p(alpha_poch + cache_size - 1, x, pol);
+           for (int i = cache_size - 1; i > 0; --i)
+           {
+              gamma_cache[i - 1] = gamma_cache[i] >= 1 ? 1 : gamma_cache[i] + regularised_gamma_prefix(alpha_poch + i - 1, x, pol, lanczos_type()) / (alpha_poch + i - 1);
+           }
+        }
         T delta_poch, alpha_poch, x, term;
+        T gamma_cache[cache_size];
         int k;
         int log_scaling;
+        int cache_offset;
         Policy pol;
      };
 
diff --git a/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp b/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp
index dba1c8eb0..ebfce2069 100644
--- a/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp
+++ b/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp
@@ -237,7 +237,6 @@
      BOOST_ASSERT(leading_a_shift > 1);
      BOOST_ASSERT(a_b_shift > 1);
 
-     boost::uintmax_t max_iter = 200;
      T first, second;
      first = boost::math::detail::hypergeometric_1F1_imp(T(a + a_shift), T(b + b_shift), z, pol, log_scaling);
      //
diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp
index 12c8a9ec6..9b31f8533 100644
--- a/test/test_1F1.hpp
+++ b/test/test_1F1.hpp
@@ -25,6 +25,10 @@
 #include <boost/math/special_functions/hypergeometric_1F1.hpp>
 #include <boost/math/quadrature/exp_sinh.hpp>
 
+#ifdef BOOST_MSVC
+#pragma warning(disable:4127)
+#endif
+
 template <class Real, class T>
 void do_test_1F1(const T& data, const char* type_name, const char* test_name)
 {
@@ -96,7 +100,7 @@ void test_spots4(T, const char* type_name)
 template <class T>
 void test_spots5(T, const char* type_name)
 {
-   std::cout << "Testing special cases" << std::endl;
+   std::cout << "Testing special cases for type " << type_name << std::endl;
    //
    // Special cases:
    //