Initial implementation of fisher_f degrees of freedom finder

include/boost/math/distributions/fisher_f.hpp

@@ -17,7 +17,84 @@
 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
 #include <boost/math/special_functions/fpclassify.hpp>
 
-namespace boost{ namespace math{
+namespace boost{ namespace math{ 
+   namespace detail{
+      template <class RealType, class Policy>
+      struct fisher_degrees_of_freedom_finder
+      {
+         fisher_degrees_of_freedom_finder(
+            RealType x_, RealType v_, bool find_v1_, RealType p_, bool c)
+            : x(x_), v(v_), find_v1(find_v1_), p(p_), comp(c) {}
+
+         RealType operator()(const RealType& input_v)
+         {
+            RealType v1 = find_v1 ? input_v : v; 
+            RealType v2 = find_v1 ? v : input_v; 
+            fisher_f_distribution<RealType, Policy> d(v1, v2);
+            return comp ?
+               p - cdf(complement(d, x))
+               : cdf(d, x) - p;
+         }
+      private:
+         RealType x;
+         RealType v;
+         bool find_v1;
+         RealType p;
+         bool comp;
+      };
+
+      template <class RealType, class Policy>
+      inline RealType find_degrees_of_freedom_fisher_f(
+         const RealType x, const RealType v, const bool find_v1, const RealType p, const RealType q, const Policy& pol)
+      {
+         BOOST_MATH_STD_USING
+         using std::fabs;
+         const char* function = "fisher_f_distribution<%1%>::find_degrees_of_freedom_f";
+         if((p == 0) || (q == 0))
+         {
+            //
+            // Can't find a thing if one of p and q is zero:
+            //
+            return policies::raise_domain_error<RealType>(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
+               RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
+         }
+         if (x < tools::epsilon<RealType>())
+         {
+            return policies::raise_evaluation_error<RealType>(function, "Can't find degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!",
+               RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
+         }
+         fisher_degrees_of_freedom_finder<RealType, Policy> f(x, v, find_v1, p < q ? p : q, p < q ? false : true);
+
+         // There are times when f has two roots - thus, two degrees of freedom can
+         // be found. We find this case by checking if the sign of f for large and 
+         // small values of v have the same sign. If the sign is the same, then there 
+         // are an even number of roots. If the signs differ, there is only one root
+         // and we can safely find the root.
+         RealType vLarge = sqrt(boost::math::tools::max_value<RealType>());
+         RealType vSmall = 1 / vLarge;
+
+         if ((f(vLarge) < 0) == (f(vSmall) < 0)){
+            return policies::raise_evaluation_error<RealType>(function, "Can't find degrees of freedom because two degrees of freedom can be found using the given parameters",
+               RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE 
+         }
+
+         tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
+         std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
+         //
+         // Pick an initial guess:
+         //
+         RealType guess = 1;
+         std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
+            f, guess, RealType(2), f(guess) < 0 ? true : false, tol, max_iter, pol);
+         RealType result = ir.first + (ir.second - ir.first) / 2;
+         if(max_iter >= policies::get_max_root_iterations<Policy>())
+         {
+            return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
+               " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE
+         }
+         return result;
+      }
+   }
 
 template <class RealType = double, class Policy = policies::policy<> >
 class fisher_f_distribution
@@ -44,7 +121,92 @@ class fisher_f_distribution
    {
       return m_df2;
    }
-
+   BOOST_MATH_GPU_ENABLED static RealType find_v1(const RealType x, const RealType v2, const RealType p)
+   {
+      constexpr auto function = "fisher_f_distribution<%1%>::find_v1";
+      typedef typename policies::evaluation<RealType, Policy>::type eval_type;
+      typedef typename policies::normalise<
+         Policy,
+         policies::promote_float<false>,
+         policies::promote_double<false>,
+         policies::discrete_quantile<>,
+         policies::assert_undefined<> >::type forwarding_policy;
+      eval_type result = detail::find_degrees_of_freedom_fisher_f(
+         static_cast<eval_type>(x),
+         static_cast<eval_type>(v2),
+         true,
+         static_cast<eval_type>(p),
+         static_cast<eval_type>(1-p),
+         forwarding_policy());
+      return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+         result,
+         function);
+   }
+   template <class A, class B, class C>
+   BOOST_MATH_GPU_ENABLED static RealType find_v1(const complemented3_type<A,B,C>& c)
+   {
+      constexpr auto function = "fisher_f_distribution<%1%>::find_v1";
+      typedef typename policies::evaluation<RealType, Policy>::type eval_type;
+      typedef typename policies::normalise<
+         Policy,
+         policies::promote_float<false>,
+         policies::promote_double<false>,
+         policies::discrete_quantile<>,
+         policies::assert_undefined<> >::type forwarding_policy;
+      eval_type result = detail::find_degrees_of_freedom_fisher_f(
+         static_cast<eval_type>(c.dist),
+         static_cast<eval_type>(c.param1),
+         true,
+         static_cast<eval_type>(1-c.param2),
+         static_cast<eval_type>(c.param2),
+         forwarding_policy());
+      return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+         result,
+         function);
+   }
+   BOOST_MATH_GPU_ENABLED static RealType find_v2(const RealType x, const RealType v2, const RealType p)
+   {
+      constexpr auto function = "fisher_f_distribution<%1%>::find_v2";
+      typedef typename policies::evaluation<RealType, Policy>::type eval_type;
+      typedef typename policies::normalise<
+         Policy,
+         policies::promote_float<false>,
+         policies::promote_double<false>,
+         policies::discrete_quantile<>,
+         policies::assert_undefined<> >::type forwarding_policy;
+      eval_type result = detail::find_degrees_of_freedom_fisher_f(
+         static_cast<eval_type>(x),
+         static_cast<eval_type>(v2),
+         false,
+         static_cast<eval_type>(p),
+         static_cast<eval_type>(1-p),
+         forwarding_policy());
+      return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+         result,
+         function);
+   }
+   template <class A, class B, class C>
+   BOOST_MATH_GPU_ENABLED static RealType find_v2(const complemented3_type<A,B,C>& c)
+   {
+      constexpr auto function = "fisher_f_distribution<%1%>::find_v2";
+      typedef typename policies::evaluation<RealType, Policy>::type eval_type;
+      typedef typename policies::normalise<
+         Policy,
+         policies::promote_float<false>,
+         policies::promote_double<false>,
+         policies::discrete_quantile<>,
+         policies::assert_undefined<> >::type forwarding_policy;
+      eval_type result = detail::find_degrees_of_freedom_fisher_f(
+         static_cast<eval_type>(c.dist),
+         static_cast<eval_type>(c.param1),
+         false,
+         static_cast<eval_type>(1-c.param2),
+         static_cast<eval_type>(c.param2),
+         forwarding_policy());
+      return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+         result,
+         function);
+   }
 private:
    //
    // Data members:

include/boost/math/distributions/non_central_f.hpp

@@ -287,7 +287,7 @@ namespace boost
          template <class A, class B, class C, class D>
          BOOST_MATH_GPU_ENABLED static RealType find_v1(const complemented4_type<A,B,C, D>& c)
          {
-            constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality";
+            constexpr auto function = "non_central_f_distribution<%1%>::find_v1";
             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
             typedef typename policies::normalise<
                Policy,
@@ -307,9 +307,9 @@ namespace boost
                result,
                function);
          }
-                  BOOST_MATH_GPU_ENABLED static RealType find_v2(const RealType x, const RealType v2, const RealType nc, const RealType p)
+         BOOST_MATH_GPU_ENABLED static RealType find_v2(const RealType x, const RealType v2, const RealType nc, const RealType p)
          {
-            constexpr auto function = "non_central_f_distribution<%1%>::find_v1";
+            constexpr auto function = "non_central_f_distribution<%1%>::find_v2";
             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
             typedef typename policies::normalise<
                Policy,
@@ -332,7 +332,7 @@ namespace boost
          template <class A, class B, class C, class D>
          BOOST_MATH_GPU_ENABLED static RealType find_v2(const complemented4_type<A,B,C, D>& c)
          {
-            constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality";
+            constexpr auto function = "non_central_f_distribution<%1%>::find_v2";
             typedef typename policies::evaluation<RealType, Policy>::type eval_type;
             typedef typename policies::normalise<
                Policy,

test/test_fisher_f.cpp

@@ -413,7 +413,19 @@ void test_spots(RealType)
     BOOST_CHECK_CLOSE(
        kurtosis(dist2)
        , static_cast<RealType>(6272) * 12 / 3456 + 3, tol2);
-    // special cases:
+
+    if (!std::is_same<RealType, boost::math::concepts::real_concept>::value)
+    {
+      RealType df1 = 10;
+      RealType df2 = 2;
+      RealType x = 6;
+      RealType P = cdf(fisher_f_distribution<RealType>(df1, df2), x);
+      BOOST_CHECK_CLOSE(fisher_f_distribution<RealType>::find_v2(x, df1, P), df2, tol2*10);
+      BOOST_CHECK_CLOSE(fisher_f_distribution<RealType>::find_v1(x, df2, P), df1, tol2*10);
+      BOOST_CHECK_CLOSE(fisher_f_distribution<RealType>::find_v1(boost::math::complement(x, df2, 1-P)), df1, tol2*10);
+      BOOST_CHECK_CLOSE(fisher_f_distribution<RealType>::find_v2(boost::math::complement(x, df1, 1-P)), df2, tol2*10);
+    }
+      // special cases:
     BOOST_MATH_CHECK_THROW(
        pdf(
           fisher_f_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),