Test beta dist for F64, F32, and F16

include/boost/math/distributions/beta.hpp

@@ -2,6 +2,7 @@
 
 // Copyright John Maddock 2006.
 // Copyright Paul A. Bristow 2006.
+// Copyright Matt Borland 2023.
 
 // Use, modification and distribution are subject to the
 // Boost Software License, Version 1.0.
@@ -241,7 +242,7 @@ namespace boost
         {
           return result;
         }
-        return ibeta_inva(beta, x, probability, Policy());
+        return static_cast<RealType>(ibeta_inva(beta, x, probability, Policy()));
       } // RealType find_alpha(beta, a, probability)
 
       static RealType find_beta(
@@ -264,7 +265,7 @@ namespace boost
         {
           return result;
         }
-        return ibeta_invb(alpha, x, probability, Policy());
+        return static_cast<RealType>(ibeta_invb(alpha, x, probability, Policy()));
       } //  RealType find_beta(alpha, x, probability)
 
     private:
@@ -396,7 +397,7 @@ namespace boost
       {
         if (a == 1)
         {
-          return 1 / beta(a, b);
+          return static_cast<RealType>(1 / beta(a, b));
         }
         else if (a < 1)
         {
@@ -411,7 +412,7 @@ namespace boost
       {
         if (b == 1)
         {
-          return 1 / beta(a, b);
+          return static_cast<RealType>(1 / beta(a, b));
         }
         else if (b < 1)
         {
@@ -423,7 +424,7 @@ namespace boost
         }
       }
       
-      return ibeta_derivative(a, b, x, Policy());
+      return static_cast<RealType>(ibeta_derivative(a, b, x, Policy()));
     } // pdf
 
     template <class RealType, class Policy>
@@ -454,7 +455,7 @@ namespace boost
       {
         return 1;
       }
-      return ibeta(a, b, x, Policy());
+      return static_cast<RealType>(ibeta(a, b, x, Policy()));
     } // beta cdf
 
     template <class RealType, class Policy>
@@ -481,16 +482,16 @@ namespace boost
       }
       if (x == 0)
       {
-        return 1;
+        return RealType(1);
       }
       else if (x == 1)
       {
-        return 0;
+        return RealType(0);
       }
       // Calculate cdf beta using the incomplete beta function.
       // Use of ibeta here prevents cancellation errors in calculating
       // 1 - x if x is very small, perhaps smaller than machine epsilon.
-      return ibetac(a, b, x, Policy());
+      return static_cast<RealType>(ibetac(a, b, x, Policy()));
     } // beta cdf
 
     template <class RealType, class Policy>
@@ -519,13 +520,13 @@ namespace boost
       // Special cases:
       if (p == 0)
       {
-        return 0;
+        return RealType(0);
       }
       if (p == 1)
       {
-        return 1;
+        return RealType(1);
       }
-      return ibeta_inv(a, b, p, static_cast<RealType*>(nullptr), Policy());
+      return static_cast<RealType>(ibeta_inv(a, b, p, static_cast<RealType*>(nullptr), Policy()));
     } // quantile
 
     template <class RealType, class Policy>
@@ -555,14 +556,14 @@ namespace boost
       // Special cases:
       if(q == 1)
       {
-        return 0;
+        return RealType(0);
       }
       if(q == 0)
       {
-        return 1;
+        return RealType(1);
       }
 
-      return ibetac_inv(a, b, q, static_cast<RealType*>(nullptr), Policy());
+      return static_cast<RealType>(ibetac_inv(a, b, q, static_cast<RealType*>(nullptr), Policy()));
     } // Quantile Complement
 
   } // namespace math

include/boost/math/special_functions/beta.hpp

@@ -1501,8 +1501,9 @@ template <class RT1, class RT2, class A>
 inline typename tools::promote_args<RT1, RT2, A>::type
    beta(RT1 a, RT2 b, A arg)
 {
-   typedef typename policies::is_policy<A>::type tag;
-   return boost::math::detail::beta(a, b, arg, static_cast<tag*>(nullptr));
+   using tag = typename policies::is_policy<A>::type;
+   using ReturnType = tools::promote_args_t<RT1, RT2, A>;
+   return static_cast<ReturnType>(boost::math::detail::beta(a, b, arg, static_cast<tag*>(nullptr)));
 }
 
 template <class RT1, class RT2>

include/boost/math/special_functions/log1p.hpp

@@ -136,24 +136,24 @@ T log1p_imp(T const& x, const Policy& pol, const std::integral_constant<int, 53>
    // Maximum Relative Change in Control Points:   8.138e-004
    // Max Error found at double precision =        3.250766e-016
    static const T P[] = {
-       0.15141069795941984e-16L,
-       0.35495104378055055e-15L,
-       0.33333333333332835L,
-       0.99249063543365859L,
-       1.1143969784156509L,
-       0.58052937949269651L,
-       0.13703234928513215L,
-       0.011294864812099712L
+       static_cast<T>(0.15141069795941984e-16L),
+       static_cast<T>(0.35495104378055055e-15L),
+       static_cast<T>(0.33333333333332835L),
+       static_cast<T>(0.99249063543365859L),
+       static_cast<T>(1.1143969784156509L),
+       static_cast<T>(0.58052937949269651L),
+       static_cast<T>(0.13703234928513215L),
+       static_cast<T>(0.011294864812099712L)
      };
    static const T Q[] = {
-       1L,
-       3.7274719063011499L,
-       5.5387948649720334L,
-       4.159201143419005L,
-       1.6423855110312755L,
-       0.31706251443180914L,
-       0.022665554431410243L,
-       -0.29252538135177773e-5L
+       static_cast<T>(1L),
+       static_cast<T>(3.7274719063011499L),
+       static_cast<T>(5.5387948649720334L),
+       static_cast<T>(4.159201143419005L),
+       static_cast<T>(1.6423855110312755L),
+       static_cast<T>(0.31706251443180914L),
+       static_cast<T>(0.022665554431410243L),
+       static_cast<T>(-0.29252538135177773e-5L)
      };
 
    T result = 1 - x / 2 + tools::evaluate_polynomial(P, x) / tools::evaluate_polynomial(Q, x);

test/test_beta_dist.cpp

@@ -52,6 +52,10 @@ using std::endl;
 #include <limits>
 using std::numeric_limits;
 
+#if __has_include(<stdfloat>)
+# include <stdfloat>
+#endif
+
 template <class RealType>
 void test_spot(
      RealType a,    // alpha a
@@ -135,6 +139,14 @@ void test_spots(RealType)
    cout << "epsilon = " << tolerance;
 
    tolerance *= 100000; // Note: NO * 100 because is fraction, NOT %.
+
+   #ifdef __STDCPP_FLOAT16_T__
+   if constexpr (std::is_same_v<RealType, std::float16_t>)
+   {
+      tolerance *= 100;
+   }
+   #endif
+
    cout  << ", Tolerance = " << tolerance * 100 << "%." << endl;
 
   // RealType teneps = boost::math::tools::epsilon<RealType>() * 10;
@@ -527,7 +539,14 @@ void test_spots(RealType)
    } // has_infinity
 
    // Error handling checks:
+   #ifdef __STDCPP_FLOAT16_T__
+   if constexpr (!std::is_same_v<std::float16_t, RealType>)
+   {
+      check_out_of_range<boost::math::beta_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.
+   }
+   #else
    check_out_of_range<boost::math::beta_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.
+   #endif
    // and range and non-finite.
 
    // Not needed??????
@@ -631,6 +650,17 @@ BOOST_AUTO_TEST_CASE( test_main )
    test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 #endif
 #endif
+
+#ifdef __STDCPP_FLOAT64_T__
+   test_spots(0.0F64);
+#endif
+#ifdef __STDCPP_FLOAT32_T__
+   test_spots(0.0F32);
+#endif
+#ifdef __STDCPP_FLOAT16_T__
+   test_spots(0.0F16);
+#endif
+
 } // BOOST_AUTO_TEST_CASE( test_main )
 
 /*