Merge pull request #1036 from boostorg/issue1034

Adjust recursion when using the Bessel function approximation with la…
boostorg · Oct 13, 2023 · 10857c0 · 10857c0
2 parents 1722ef9 + 57b09f4
commit 10857c0
Show file tree

Hide file tree

Showing 5 changed files with 35 additions and 7 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
@@ -393,7 +393,15 @@
         // So shift b to match a (b shifting seems to be more stable via method of ratios).
         //
         int b_shift = itrunc(b - a);
+        if ((b_shift < 0) && (b - b_shift != a))
+           b_shift -= 1;
         T b_local = b - b_shift;
+        if ((b_local - a - 0.5 <= 0) && (b_local != a))
+        {
+           // Make sure b_local - a - 0.5 > 0
+           b_shift -= 1;
+           b_local += 1;
+        }
         T h = boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b_local, z, T(b_local - a), pol, log_scaling);
         return hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling);
      }

diff --git a/include/boost/math/special_functions/hypergeometric_1F1.hpp b/include/boost/math/special_functions/hypergeometric_1F1.hpp
@@ -355,7 +355,9 @@ namespace boost { namespace math { namespace detail {
          if ((a < 0) && (a == ceil(a)) && (a > -50))
             return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);
 
-         return (b + z) * exp(z) / b;
+         log_scaling = lltrunc(floor(z));
+         T local_z = z - log_scaling;
+         return (b + z) * exp(local_z) / b;
       }
 
       if ((a == 1) && (b == 2))

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>, 183> hypergeometric_1F1_bugs = { {
+   static const std::array<std::array<T, 4>, 186> 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) }} ,
@@ -390,6 +390,11 @@ void test_spots6(T, const char* type_name)
       { { -28, 28, 28, SC_(-6.2125286411657869483728921158018766e-9) } }, 
       { { -29, 29, 29, SC_(-1.3521578972057573423569167878458172e-9) } }, 
       {{ -30, 30, 30, SC_(8.2238878884841599031462003461115991e-10) } },
+
+      // https://github.com/boostorg/math/issues/1034
+      {{ 13, 1.5f, 61, SC_(1.35508577094765660270265300640877455638098585524020525369044e39)}},
+      {{ 13, 1.5f - T(1)/128, 61, SC_(1.40067238333701986992154961431485209677766220602448290643906e39)}},
+      {{ 13, 1.5f + T(1)/128, 61, SC_(1.31105748771677778012064837998217769289913724450105998963999e39)}},
    } };
    static const std::array<std::array<T, 4>, 2> hypergeometric_1F1_big_bugs = { {
 #if DBL_MAX_EXP == LDBL_MAX_EXP

diff --git a/test/test_1F1_log.cpp b/test/test_1F1_log.cpp
@@ -160,11 +160,6 @@ BOOST_AUTO_TEST_CASE( test_main )
    expected_results();
    BOOST_MATH_CONTROL_FP;
 
-#if !defined(TEST) || (TEST == 1)
-   test_hypergeometric_mellin_transform<double>();
-   test_hypergeometric_laplace_transform<double>();
-#endif
-
 #ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS
 #if !defined(TEST) || (TEST == 2)
    test_spots(0.0F, "float");

diff --git a/test/test_1F1_log.hpp b/test/test_1F1_log.hpp
@@ -77,10 +77,28 @@ void test_spots2(T, const char* type_name)
    do_test_1F1<T>(hypergeometric_1f1_log_large_unsolved, type_name, "Large random values - log - unsolved");
 }
 
+template <class T>
+void test_spots_bugs(T, const char* type_name)
+{
+   static const std::array<std::array<T, 4>, 7> hypergeometric_1F1_bugs = { {
+    // Found while investigating https://github.com/boostorg/math/issues/1034
+    {{ 21156.0f, 21156.0f, 11322, SC_(11322.0)}},
+    {{ 21156.0f, 21155.0f, 11322, SC_(11322.428655862323560632951631114666466652986288119296800531328684)}},
+    {{ 21156.0f, 21154.0f, 11322, SC_(11322.857338938931770780542471014439235046236959098048505808665509)}},
+    {{ 21156.0f, 21154.5f, 11322, SC_(11322.642993998700652342915766513423502332460377941025163971463001)}},
+    {{ 21156.0f, 21155.0f - 1.0f / 128, 11322, SC_(11322.43200484338133063401686149227756707)}},
+    {{ 21156.0f, 21154.5f - 1.0f / 128, 11322, SC_(11322.646343086066466097278446364687150256282052775170519455915995)}},
+    {{ 21156.0f, 21154.0f + 1.0f / 128, 11322, SC_(11322.85398974691465958225700429672975704)}},
+   } };
+
+   do_test_1F1<T>(hypergeometric_1F1_bugs, type_name, "Large random values - log - bug cases");
+}
+
 template <class T>
 void test_spots(T z, const char* type_name)
 {
    test_spots1(z, type_name);
+   test_spots_bugs(z, type_name);
 #ifdef TEST_UNSOLVED
    test_spots2(z, type_name);
 #endif