diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp
index 02a75a28d..a6df68b8f 100644
--- a/include/boost/math/distributions/non_central_f.hpp
+++ b/include/boost/math/distributions/non_central_f.hpp
@@ -46,8 +46,7 @@ namespace boost
          };
          
          template <class RealType, class Policy>
-         BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const Policy& pol)
-         {
+         BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const RealType p_q_precision, const Policy& pol)         {
             constexpr auto function = "non_central_f<%1%>::find_non_centrality";
 
             if ( p == 0 || q == 0) {
@@ -56,17 +55,20 @@ namespace boost
             }
 
             // Check if nc = 0 (which is just the F-distribution)
-            // I first tried comparing to boost::math::tools::epsilon<RealType>(),
-            // but this was never satisfied for floats or doubles. Only 
-            // long doubles would correctly return 0.
-            fisher_f_distribution<RealType, Policy> dist(v1, v2); 
-            if (boost::math::relative_difference(pdf(dist, x), p) < 1e-7){
+            non_centrality_finder_f<RealType, Policy> f(x, v1, v2, p < q ? p : q, p < q ? false : true);
+            // This occurs when the root finder would need to find a result smaller than
+            // tools::min_value (which it cannot do).  Note that we have to add in a small
+            // amount of "tolerance" since the subtraction in our termination condition
+            // implies a small amount of wobble in the result which should be of the
+            // order p * eps (note not q * eps, since q is calculated as 1-p).
+            // Also note that p_q_precision is passed down from our caller as the
+            // epsilon of the original called values, and not after possible promotion.
+            if (f(tools::min_value<RealType>()) <= 3 * p_q_precision * p){
                return 0;
             }
 
-            non_centrality_finder_f<RealType, Policy> f(x, v1, v2, p < q ? p : q, p < q ? false : true);
             RealType guess = RealType(10);                       // Starting guess.
-            RealType factor = RealType(2);                                 // How big steps to take when searching.
+            RealType factor = RealType(2);                       // How big steps to take when searching.
             boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
             tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
 
@@ -134,6 +136,7 @@ namespace boost
                static_cast<eval_type>(v2),
                static_cast<eval_type>(p),
                static_cast<eval_type>(1-p),
+               static_cast<eval_type>(tools::epsilon<RealType>()),
                forwarding_policy());
             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                result,
@@ -156,6 +159,7 @@ namespace boost
                static_cast<eval_type>(c.param2),
                static_cast<eval_type>(1-c.param3),
                static_cast<eval_type>(c.param3),
+               static_cast<eval_type>(tools::epsilon<RealType>()),
                forwarding_policy());
             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                result,
diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp
index 3ccb9e3b4..6e8da36e9 100644
--- a/test/test_nc_f.cpp
+++ b/test/test_nc_f.cpp
@@ -159,12 +159,10 @@ void test_spot(
 }
 
 template <class RealType> // Any floating-point type RealType.
-void test_spots(RealType)
+void test_spots(RealType, const char* name = nullptr)
 {
    RealType tolerance = boost::math::tools::epsilon<RealType>() * 10000;
-
-   cout << "Tolerance = " << (tolerance / 100) << "%." << endl;
-
+   std::cout << "Testing spot values with type " << name << " (Tolerance = " << (tolerance / 100) << "%)." << std::endl;
    //
    // Spot tests from Mathematica computed values:
    //
@@ -337,8 +335,8 @@ void test_spots(RealType)
    boost::math::non_central_f_distribution<RealType> dist_no_centrality(a, b, nc);
    for (RealType x : x_vals)
    {
-      RealType P = pdf(dist_no_centrality, x);
-      BOOST_CHECK_CLOSE(dist.find_non_centrality(x, a, b, P), static_cast<RealType>(0), tolerance);
+      RealType P = cdf(dist_no_centrality, x);
+      BOOST_CHECK(dist.find_non_centrality(x, a, b, P) < tolerance);
    }
    // Case when P=1 or P=0 
    BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 1), std::domain_error);
@@ -354,12 +352,12 @@ BOOST_AUTO_TEST_CASE( test_main )
    // Basic sanity-check spot values.
    expected_results();
    // (Parameter value, arbitrarily zero, only communicates the floating point type).
-   test_spots(0.0F); // Test float.
-   test_spots(0.0); // Test double.
+   test_spots(0.0F, "float"); // Test float.
+   test_spots(0.0, "double"); // Test double.
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-   test_spots(0.0L); // Test long double.
+   test_spots(0.0L, "long double"); // Test long double.
 #if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS)
-   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
+   test_spots(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
 #endif
 #endif