Merge pull request #778 from boostorg/hyper1f1_bugfix

Minor 1F1 bug fixes.
boostorg · Mar 13, 2022 · 2d8413e · 2d8413e
2 parents 0bb12b2 + bd7fcf1
commit 2d8413e
Show file tree

Hide file tree

Showing 3 changed files with 28 additions and 4 deletions.
diff --git a/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp b/include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp
@@ -291,6 +291,14 @@
        ak += 2;
        integer_part -= 2;
     }
+    if (ak - 1 == b)
+    {
+       // When ak - 1 == b are recursion coefficients dissappear to zero and
+       // we end up with a NaN result.  Reduce the recursion steps by 1 to
+       // avoid this.  We rely on |b| small and therefore no infinite recursion.
+       ak -= 1;
+       integer_part += 1;
+    }
 
     if (-integer_part > static_cast<std::intmax_t>(policies::get_max_series_iterations<Policy>()))
        return policies::raise_evaluation_error<T>(function, "1F1 arguments sit in a range with a so negative that we have no evaluation method, got a = %1%", std::numeric_limits<T>::quiet_NaN(), pol);

diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp
@@ -132,6 +132,7 @@
         Real scaling_factor = exp(Real(log_scaling_factor));
         Real term_m1;
         long long local_scaling = 0;
+        bool have_no_correct_bits = false;
 
         if ((aj.size() == 1) && (bj.size() == 0))
         {
@@ -238,10 +239,20 @@
               break;
            if (abs_result * tol > abs(result))
            {
-              // We have no correct bits in the result... just give up!
-              result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol);
-              return std::make_pair(result, result);
+              // Check if result is so small compared to abs_resuslt that there are no longer any
+              // correct bits... we require two consecutive passes here before aborting to
+              // avoid false positives when result transiently drops to near zero then rebounds.
+              if (have_no_correct_bits)
+              {
+                 // We have no correct bits in the result... just give up!
+                 result = boost::math::policies::raise_evaluation_error("boost::math::hypergeometric_pFq<%1%>", "Cancellation is so severe that no bits in the reuslt are correct, last result was %1%", Real(result * exp(Real(log_scale))), pol);
+                 return std::make_pair(result, result);
+              }
+              else
+                 have_no_correct_bits = true;
            }
+           else
+              have_no_correct_bits = false;
            term0 = term;
         }
         //std::cout << "result = " << result << std::endl;

diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp
@@ -162,7 +162,7 @@ void test_spots5(T, const char* type_name)
 template <class T>
 void test_spots6(T, const char* type_name)
 {
-   static const std::array<std::array<T, 4>, 91> hypergeometric_1F1_bugs = { {
+   static const std::array<std::array<T, 4>, 93> hypergeometric_1F1_bugs = { {
         { { static_cast<double>(17955.561660766602), static_cast<double>(9.6968994205831605e-09), static_cast<double>(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }},
         { { static_cast<double>(17955.561660766602), static_cast<double>(-9.6968994205831605e-09), static_cast<double>(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }},
         { { static_cast<double>(-17955.561660766602), static_cast<double>(-9.6968994205831605e-09), static_cast<double>(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} ,
@@ -287,6 +287,11 @@ void test_spots6(T, const char* type_name)
         { { (T)std::ldexp((double)-9751199809536000, -45), (T)std::ldexp((double)-17654191685632000, -47), (T)std::ldexp((double)10587451850752000, -47), SC_(-1.9601415510439595625538337964298353914980331018955e+68) }},
         { { (T)std::ldexp((double)-15233620754432000, -45), (T)std::ldexp((double)-12708283072512000, -46), (T)std::ldexp((double)10255461007360000, -46), SC_(-5.4344106361679075861859567858016187271235441673635e+125) }},
         { { (T)std::ldexp((double)-11241354149888000, -45), (T)std::ldexp((double)-9580579905536000, -45), (T)std::ldexp((double)12224976846848000, -47), SC_(12046856548470067405870726490464935201150430438.035) }},
+        //
+        // Bugs found while testing color maps:
+        //
+        { { SC_(0.078125000000000000), SC_(-0.039062500000000000), SC_(0.5), SC_(-0.3371910410915676603577770246237158427221) }},
+        { { SC_(-19.492187500000000), SC_(0.50781250000000000), SC_(0.5), SC_(1.2551298228307647570646714060395253358015) }},
    } };
    static const std::array<std::array<T, 4>, 2> hypergeometric_1F1_big_bugs = { {
 #if DBL_MAX_EXP == LDBL_MAX_EXP