Improve erf/expm1/expint coverage.

include/boost/math/special_functions/detail/bessel_ik.hpp

@@ -326,6 +326,9 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
     T scale = 1;
     T scale_sign = 1;
 
+    n = iround(v, pol);
+    u = v - n;                              // -1/2 <= u < 1/2
+
     if (((kind & need_i) == 0) && (fabs(4 * v * v - 25) / (8 * x) < tools::forth_root_epsilon<T>()))
     {
        // A&S 9.7.2
@@ -337,8 +340,6 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
     }
     else
     {
-       n = iround(v, pol);
-       u = v - n;                              // -1/2 <= u < 1/2
        BOOST_MATH_INSTRUMENT_VARIABLE(n);
        BOOST_MATH_INSTRUMENT_VARIABLE(u);
 
@@ -412,7 +413,7 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol)
        else
           Iv = std::numeric_limits<T>::quiet_NaN(); // any value will do
     }
-    if (reflect)
+    if (reflect && (kind & need_i))
     {
         T z = (u + n % 2);
         T fact = (2 / pi<T>()) * (boost::math::sin_pi(z, pol) * Kv);

include/boost/math/special_functions/detail/bessel_j1.hpp

@@ -34,36 +34,9 @@ namespace boost { namespace math{  namespace detail{
 template <typename T>
 T bessel_j1(T x);
 
-template <class T>
-struct bessel_j1_initializer
-{
-   struct init
-   {
-      init()
-      {
-         do_init();
-      }
-      static void do_init()
-      {
-         bessel_j1(T(1));
-      }
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T>
-const typename bessel_j1_initializer<T>::init bessel_j1_initializer<T>::initializer;
-
 template <typename T>
 T bessel_j1(T x)
 {
-    bessel_j1_initializer<T>::force_instantiate();
-
     static const T P1[] = {
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)),
          static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)),

include/boost/math/special_functions/detail/bessel_k0.hpp

@@ -46,40 +46,6 @@ namespace boost { namespace math { namespace detail{
 template <typename T>
 T bessel_k0(const T& x);
 
-template <class T, class tag>
-struct bessel_k0_initializer
-{
-   struct init
-   {
-      init()
-      {
-         do_init(tag());
-      }
-      static void do_init(const std::integral_constant<int, 113>&)
-      {
-         bessel_k0(T(0.5));
-         bessel_k0(T(1.5));
-      }
-      static void do_init(const std::integral_constant<int, 64>&)
-      {
-         bessel_k0(T(0.5));
-         bessel_k0(T(1.5));
-      }
-      template <class U>
-      static void do_init(const U&){}
-      void force_instantiate()const{}
-   };
-   static const init initializer;
-   static void force_instantiate()
-   {
-      initializer.force_instantiate();
-   }
-};
-
-template <class T, class tag>
-const typename bessel_k0_initializer<T, tag>::init bessel_k0_initializer<T, tag>::initializer;
-
-
 template <typename T, int N>
 T bessel_k0_imp(const T&, const std::integral_constant<int, N>&)
 {
@@ -505,7 +471,6 @@ inline T bessel_k0(const T& x)
       113 : -1
    > tag_type;
 
-   bessel_k0_initializer<T, tag_type>::force_instantiate();
    return bessel_k0_imp(x, tag_type());
 }
 

include/boost/math/special_functions/detail/bessel_k1.hpp

@@ -46,41 +46,6 @@ namespace boost { namespace math { namespace detail{
    template <typename T>
    T bessel_k1(const T&);
 
-   template <class T, class tag>
-   struct bessel_k1_initializer
-   {
-      struct init
-      {
-         init()
-         {
-            do_init(tag());
-         }
-         static void do_init(const std::integral_constant<int, 113>&)
-         {
-            bessel_k1(T(0.5));
-            bessel_k1(T(2));
-            bessel_k1(T(6));
-         }
-         static void do_init(const std::integral_constant<int, 64>&)
-         {
-            bessel_k1(T(0.5));
-            bessel_k1(T(6));
-         }
-         template <class U>
-         static void do_init(const U&) {}
-         void force_instantiate()const {}
-      };
-      static const init initializer;
-      static void force_instantiate()
-      {
-         initializer.force_instantiate();
-      }
-   };
-
-   template <class T, class tag>
-   const typename bessel_k1_initializer<T, tag>::init bessel_k1_initializer<T, tag>::initializer;
-
-
    template <typename T, int N>
    inline T bessel_k1_imp(const T&, const std::integral_constant<int, N>&)
    {
@@ -547,7 +512,6 @@ namespace boost { namespace math { namespace detail{
          113 : -1
       > tag_type;
 
-      bessel_k1_initializer<T, tag_type>::force_instantiate();
       return bessel_k1_imp(x, tag_type());
    }
 

include/boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp

@@ -20,12 +20,11 @@
   {
     BOOST_MATH_STD_USING
 
-    const bool is_z_nonpositive = z <= 0;
+    //const bool is_z_nonpositive = z <= 0;
+    BOOST_MATH_ASSERT(z < 0);  // condition used at call site
 
-    const T sqrt_z = is_z_nonpositive ? T(sqrt(-z)) : T(sqrt(z));
-    const T bessel_mult = is_z_nonpositive ?
-      boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol) :
-      boost::math::cyl_bessel_i(b - 1, 2 * sqrt_z, pol) ;
+    const T sqrt_z = sqrt(-z);
+    const T bessel_mult = boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol);
 
     if (b > boost::math::max_factorial<T>::value)
     {

include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp

@@ -74,8 +74,12 @@
            {
               // We get very limited precision due to rapid denormalisation/underflow of the Bessel values, raise an exception and try something else:
               policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Underflow in Bessel functions", bessel_cache[cache_size - 1], pol);
+              // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later:
+              std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits<T>::quiet_NaN());
+              cache_offset = -cache_size;
+              return;
            }
-           if ((term * bessel_cache[cache_size - 1] < tools::min_value<T>() / (tools::epsilon<T>() * tools::epsilon<T>())) || !(boost::math::isfinite)(term) || (!std::numeric_limits<T>::has_infinity && (fabs(term) > tools::max_value<T>())))
+           if ((fabs(term * bessel_cache[cache_size - 1]) < tools::min_value<T>() / (tools::epsilon<T>() * tools::epsilon<T>())) || !(boost::math::isfinite)(term) || (!std::numeric_limits<T>::has_infinity && (fabs(term) > tools::max_value<T>())))
            {
               term = -log(fabs(bessel_arg)) * b_minus_1_plus_n / 2;
               log_scale = lltrunc(term);
@@ -88,15 +92,27 @@
            if constexpr (std::numeric_limits<T>::has_infinity)
            {
               if (!(boost::math::isfinite)(bessel_cache[cache_size - 1]))
+              {
                  policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol);
+                 // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later:
+                 std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits<T>::quiet_NaN());
+              }
            }
            else
               if ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))
+              {
                  policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol);
+                 // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later:
+                 std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits<T>::quiet_NaN());
+              }
 #else
            if ((std::numeric_limits<T>::has_infinity && !(boost::math::isfinite)(bessel_cache[cache_size - 1])) 
               || (!std::numeric_limits<T>::has_infinity && ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))))
+           {
               policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol);
+              // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later:
+              std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits<T>::quiet_NaN());
+           }
 #endif
            cache_offset = -cache_size;
            refill_cache();
@@ -108,8 +124,13 @@
            // very small (or zero) when b == 2a:
            //
            BOOST_MATH_STD_USING
+           //
+           // Except in the multiprecision case, we have probably illiminated anything
+           // would need more than the default 64 Bessel Functions.  Anything more
+           // than that risks becoming a divergent series anyway...
+           //
            if(n - 2 - cache_offset >= cache_size)
-              refill_cache();
+              refill_cache();  // LCOV_EXCL_LINE
            T result = A_minus_2 * term * bessel_cache[n - 2 - cache_offset];
            term *= mult;
            ++n;
@@ -122,7 +143,7 @@
            if (A_minus_2 != 0)
            {
               if (n - 2 - cache_offset >= cache_size)
-                 refill_cache();
+                 refill_cache(); // LCOV_EXCL_LINE
               result += A_minus_2 * term * bessel_cache[n - 2 - cache_offset];
            }
            term *= mult;

include/boost/math/special_functions/detail/igamma_large.hpp

@@ -88,6 +88,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 64
 
    T workspace[13];
 
+   // LCOV_EXCL_START
    static const T C0[] = {
       BOOST_MATH_BIG_CONSTANT(T, 64, -0.333333333333333333333),
       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0833333333333333333333),
@@ -265,6 +266,8 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 64
       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00640336283380806979482),
       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00404101610816766177474),
    };
+   // LCOV_EXCL_END
+
    workspace[12] = tools::evaluate_polynomial(C12, z);
 
    T result = tools::evaluate_polynomial<13, T, T>(workspace, 1/a);
@@ -293,6 +296,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 53
 
    T workspace[10];
 
+   // LCOV_EXCL_START
    static const T C0[] = {
       static_cast<T>(-0.33333333333333333L),
       static_cast<T>(0.083333333333333333L),
@@ -406,6 +410,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 53
       static_cast<T>(0.00083949872067208728L),
       static_cast<T>(-0.00043829709854172101L),
    };
+   // LCOV_EXCL_END
    workspace[8] = tools::evaluate_polynomial(C8, z);
    workspace[9] = static_cast<T>(-0.00059676129019274625L);
 
@@ -435,6 +440,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 24
 
    T workspace[3];
 
+   // LCOV_EXCL_START
    static const T C0[] = {
       static_cast<T>(-0.333333333L),
       static_cast<T>(0.0833333333L),
@@ -461,6 +467,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 24
       static_cast<T>(0.000771604938L),
    };
    workspace[2] = tools::evaluate_polynomial(C2, z);
+   // LCOV_EXCL_END
 
    T result = tools::evaluate_polynomial(workspace, 1/a);
    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
@@ -478,6 +485,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 24
 // It's use for a < 200 is not recommended, that would
 // require many more terms in the polynomials.
 //
+// LCOV_EXCL_START: 128-bit floats not deliberately tested in our coverage tests (takes too long)
 template <class T, class Policy>
 T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 113> const *)
 {
@@ -768,6 +776,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant<int, 11
 
    return result;
 }
+// LCOV_EXCL_END
 
 }  // namespace detail
 }  // namespace math