Non Central T: track cancellation errors,

include/boost/math/distributions/non_central_t.hpp

@@ -249,8 +249,7 @@ namespace boost
                // t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31.
                // "Continuous Univariate Distributions".  N.L. Johnson, S. Kotz and
                // N. Balkrishnan. 1995. John Wiley and Sons New York.
-               T result = cdf(students_t_distribution<T, Policy>(v), t - delta);
-               return invert ? 1 - result : result;
+               return invert ? cdf(complement(students_t_distribution<T, Policy>(v), t - delta)) : cdf(students_t_distribution<T, Policy>(v), t - delta);
             }
             //
             // x and y are the corresponding random
@@ -274,17 +273,31 @@ namespace boost
                //
                // Calculate p:
                //
+               T c = 0;
                if(x != 0)
                {
-                  result = non_central_beta_p(a, b, d2, x, y, pol);
-                  result = non_central_t2_p(v, delta, x, y, pol, result);
+                  c = non_central_beta_p(a, b, d2, x, y, pol);
+                  result = non_central_t2_p(v, delta, x, y, pol, c);
                   result /= 2;
                }
                else
                   result = 0;
                if (invert)
                {
                   result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta)) - result;
+                  if ((x != 0) && (fabs(result / (c * tools::epsilon<T>())) < 1000))
+                  {
+                      // We've cancelled out most of the digits in the result, try A&S 26.7.9,
+                      // this is only accurate to a couple of digits at best, but it's better
+                      // than nothing when all else fails, see https://github.com/boostorg/math/issues/1430
+                      t = -t;
+                      delta = -delta;
+
+                      T z = (t * (1 - 1 / (4 * v)) - delta) / sqrt(1 + t * t / (2 * v));
+
+                      return cdf(boost::math::normal_distribution<T, Policy>(), z);
+
+                  }
                   invert = false;
                }
                else

test/test_nc_t.hpp

@@ -344,6 +344,11 @@ void test_spots(RealType)
       distro1 d(8.0f, 8.5f);
       BOOST_CHECK_CLOSE(pdf(d, -1), static_cast<RealType>(6.1747948083757028903541988987716621647020752431287e-20), 2e-5);  // Can we do better on accuracy here?
    }
+   // https://github.com/boostorg/math/issues/1430
+   {
+       distro1 d(1000.f, 23.f);
+       BOOST_CHECK_CLOSE_FRACTION(cdf(d, -1), static_cast<RealType>(1.61471461239552e-127), 1e-3);
+   }
 
 } // template <class RealType>void test_spots(RealType)
 
@@ -542,8 +547,10 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test)
             BOOST_ERROR(e.what());
          }
          // Code coverage:
-         BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], boost::math::tools::epsilon<value_type>() / 2, data[i][3]), boost::math::evaluation_error);
-         BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], boost::math::tools::epsilon<value_type>() / 2, data[i][3])), boost::math::evaluation_error);
+         using no_promote_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false> >;
+         using no_promote_distro = boost::math::non_central_t_distribution<value_type, no_promote_policy>;
+         BOOST_CHECK_THROW(no_promote_distro::find_degrees_of_freedom(data[i][1], boost::math::tools::epsilon<value_type>() / 2, data[i][3]), boost::math::evaluation_error);
+         BOOST_CHECK_THROW(no_promote_distro::find_degrees_of_freedom(boost::math::complement(data[i][1], boost::math::tools::epsilon<value_type>() / 2, data[i][3])), boost::math::evaluation_error);
          BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], value_type(0)), boost::math::evaluation_error);
          BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], value_type(0))), boost::math::evaluation_error);
          BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], value_type(1)), boost::math::evaluation_error);