Add manual promotion to arcsine, beta, and Bernoulli distributions

include/boost/math/distributions/arcsine.hpp

@@ -239,7 +239,12 @@ namespace boost
       {
         return result;
       }
-      return  (x_min + x_max) / 2;
+
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      const auto promoted_x_max = static_cast<policy_promoted_type>(x_max);
+      const auto promoted_x_min = static_cast<policy_promoted_type>(x_min);
+
+      return static_cast<RealType>((promoted_x_min + promoted_x_max) / 2);
     } // mean
 
     template <class RealType, class Policy>
@@ -257,7 +262,12 @@ namespace boost
       {
         return result;
       }
-      return  (x_max - x_min) * (x_max - x_min) / 8;
+
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      const auto promoted_x_max = static_cast<policy_promoted_type>(x_max);
+      const auto promoted_x_min = static_cast<policy_promoted_type>(x_min);
+
+      return static_cast<RealType>((promoted_x_max - promoted_x_min) * (promoted_x_max - promoted_x_min) / 8);
     } // variance
 
     template <class RealType, class Policy>
@@ -286,7 +296,12 @@ namespace boost
       {
         return result;
       }
-      return  (x_min + x_max) / 2;
+
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      const auto promoted_x_max = static_cast<policy_promoted_type>(x_max);
+      const auto promoted_x_min = static_cast<policy_promoted_type>(x_min);
+
+      return static_cast<RealType>((promoted_x_min + promoted_x_max) / 2);
     }
 
     template <class RealType, class Policy>
@@ -324,8 +339,10 @@ namespace boost
       {
         return result;
       }
-      result = -3;
-      return  result / 2;
+
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      result = static_cast<RealType>(static_cast<policy_promoted_type>(-3) / static_cast<policy_promoted_type>(2));
+      return result;
     } // kurtosis_excess
 
     template <class RealType, class Policy>
@@ -345,7 +362,9 @@ namespace boost
         return result;
       }
 
-      return 3 + kurtosis_excess(dist);
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      result = static_cast<RealType>(3 + static_cast<policy_promoted_type>(-3) / static_cast<policy_promoted_type>(2));
+      return result;
     } // kurtosis
 
     template <class RealType, class Policy>
@@ -369,8 +388,15 @@ namespace boost
       {
         return result;
       }
+
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
       using boost::math::constants::pi;
-      result = static_cast<RealType>(1) / (pi<RealType>() * sqrt((x - lo) * (hi - x)));
+
+      const auto promoted_lo = static_cast<policy_promoted_type>(lo);
+      const auto promoted_hi = static_cast<policy_promoted_type>(hi);
+      const auto promoted_x = static_cast<policy_promoted_type>(x);
+
+      result = static_cast<RealType>(static_cast<policy_promoted_type>(1) / (pi<policy_promoted_type>() * sqrt((promoted_x - promoted_lo) * (promoted_hi - promoted_x))));
       return result;
     } // pdf
 
@@ -402,8 +428,15 @@ namespace boost
       {
         return 1;
       }
+
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
       using boost::math::constants::pi;
-      result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>();
+
+      const auto promoted_x = static_cast<policy_promoted_type>(x);
+      const auto promoted_x_max = static_cast<policy_promoted_type>(x_max);
+      const auto promoted_x_min = static_cast<policy_promoted_type>(x_min);
+
+      result = static_cast<RealType>(static_cast<policy_promoted_type>(2) * asin(sqrt((promoted_x - promoted_x_min) / (promoted_x_max - promoted_x_min))) / pi<policy_promoted_type>());
       return result;
     } // arcsine cdf
 
@@ -435,11 +468,19 @@ namespace boost
       {
         return 1;
       }
-      using boost::math::constants::pi;
+
       // Naive version x = 1 - x;
       // result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>();
       // is less accurate, so use acos instead of asin for complement.
-      result = static_cast<RealType>(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>();
+
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      using boost::math::constants::pi;
+
+      const auto promoted_x = static_cast<policy_promoted_type>(x);
+      const auto promoted_x_max = static_cast<policy_promoted_type>(x_max);
+      const auto promoted_x_min = static_cast<policy_promoted_type>(x_min);
+
+      result = static_cast<RealType>(static_cast<policy_promoted_type>(2) * acos(sqrt((promoted_x - promoted_x_min) / (promoted_x_max - promoted_x_min))) / pi<policy_promoted_type>());
       return result;
     } // arcsine ccdf
 
@@ -480,9 +521,13 @@ namespace boost
         return 1;
       }
 
-      RealType sin2hpip = sin(half_pi<RealType>() * p);
-      RealType sin2hpip2 = sin2hpip * sin2hpip;
-      result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2;
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      const policy_promoted_type sin2hpip = sin(half_pi<policy_promoted_type>() * static_cast<policy_promoted_type>(p));
+      const policy_promoted_type sin2hpip2 = sin2hpip * sin2hpip;
+      const auto promoted_x_min = static_cast<policy_promoted_type>(x_min);
+      const auto promoted_x_max = static_cast<policy_promoted_type>(x_max);
+
+      result = static_cast<RealType>(-promoted_x_min * sin2hpip2 + promoted_x_min + promoted_x_max * sin2hpip2);
 
       return result;
     } // quantile
@@ -526,9 +571,13 @@ namespace boost
       //result = cos(half_pi<RealType>() * q); // for arcsine(0,1)
       //result = result * result;
       // For generalized arcsine:
-      RealType cos2hpip = cos(half_pi<RealType>() * q);
-      RealType cos2hpip2 = cos2hpip * cos2hpip;
-      result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2;
+      using policy_promoted_type = policies::evaluation_t<RealType, Policy>;
+      const policy_promoted_type cos2hpip = cos(half_pi<policy_promoted_type>() * static_cast<policy_promoted_type>(q));
+      const policy_promoted_type cos2hpip2 = cos2hpip * cos2hpip;
+      const auto promoted_x_min = static_cast<policy_promoted_type>(x_min);
+      const auto promoted_x_max = static_cast<policy_promoted_type>(x_max);
+
+      result = -promoted_x_min * cos2hpip2 + promoted_x_min + promoted_x_max * cos2hpip2;
 
       return result;
     } // Quantile Complement

include/boost/math/distributions/bernoulli.hpp

@@ -311,26 +311,31 @@ namespace boost
     BOOST_MATH_GPU_ENABLED inline RealType skewness(const bernoulli_distribution<RealType, Policy>& dist)
     {
       BOOST_MATH_STD_USING; // Aid ADL for sqrt.
-      RealType p = dist.success_fraction();
-      return (1 - 2 * p) / sqrt(p * (1 - p));
+
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      const promoted_real_type p = static_cast<promoted_real_type>(dist.success_fraction());
+      return static_cast<RealType>((1 - 2 * p) / sqrt(p * (1 - p)));
     }
 
     template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const bernoulli_distribution<RealType, Policy>& dist)
     {
-      RealType p = dist.success_fraction();
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      promoted_real_type p = static_cast<promoted_real_type>(dist.success_fraction());
       // Note Wolfram says this is kurtosis in text, but gamma2 is the kurtosis excess,
       // and Wikipedia also says this is the kurtosis excess formula.
       // return (6 * p * p - 6 * p + 1) / (p * (1 - p));
       // But Wolfram kurtosis article gives this simpler formula for kurtosis excess:
-      return 1 / (1 - p) + 1/p -6;
+      return static_cast<RealType>(1 / (1 - p) + 1/p - 6);
     }
 
     template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const bernoulli_distribution<RealType, Policy>& dist)
     {
-      RealType p = dist.success_fraction();
-      return 1 / (1 - p) + 1/p -6 + 3;
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      promoted_real_type p = static_cast<promoted_real_type>(dist.success_fraction());
+
+      return static_cast<RealType>(1 / (1 - p) + 1/p - 6 + 3);
       // Simpler than:
       // return (6 * p * p - 6 * p + 1) / (p * (1 - p)) + 3;
     }

include/boost/math/distributions/beta.hpp

@@ -298,15 +298,21 @@ namespace boost
     template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType mean(const beta_distribution<RealType, Policy>& dist)
     { // Mean of beta distribution = np.
-      return  dist.alpha() / (dist.alpha() + dist.beta());
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      const auto a = static_cast<promoted_real_type>(dist.alpha());
+      const auto b = static_cast<promoted_real_type>(dist.beta());
+
+      return static_cast<RealType>(a / (a + b));
     } // mean
 
     template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType variance(const beta_distribution<RealType, Policy>& dist)
     { // Variance of beta distribution = np(1-p).
-      RealType a = dist.alpha();
-      RealType b = dist.beta();
-      return  (a * b) / ((a + b ) * (a + b) * (a + b + 1));
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      const auto a = static_cast<promoted_real_type>(dist.alpha());
+      const auto b = static_cast<promoted_real_type>(dist.beta());
+
+      return static_cast<RealType>((a * b) / ((a + b ) * (a + b) * (a + b + 1)));
     } // variance
 
     template <class RealType, class Policy>
@@ -330,9 +336,11 @@ namespace boost
           "mode undefined for beta = %1%, must be > 1!", dist.beta(), Policy());
         return result;
       }
-      RealType a = dist.alpha();
-      RealType b = dist.beta();
-      return (a-1) / (a + b - 2);
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      const auto a = static_cast<promoted_real_type>(dist.alpha());
+      const auto b = static_cast<promoted_real_type>(dist.beta());
+
+      return static_cast<RealType>((a-1) / (a + b - 2));
     } // mode
 
     //template <class RealType, class Policy>
@@ -347,20 +355,23 @@ namespace boost
     BOOST_MATH_GPU_ENABLED inline RealType skewness(const beta_distribution<RealType, Policy>& dist)
     {
       BOOST_MATH_STD_USING // ADL of std functions.
-      RealType a = dist.alpha();
-      RealType b = dist.beta();
-      return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b));
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      const auto a = static_cast<promoted_real_type>(dist.alpha());
+      const auto b = static_cast<promoted_real_type>(dist.beta());
+
+      return static_cast<RealType>((2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b)));
     } // skewness
 
     template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const beta_distribution<RealType, Policy>& dist)
     {
-      RealType a = dist.alpha();
-      RealType b = dist.beta();
-      RealType a_2 = a * a;
-      RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2));
-      RealType d = a * b * (a + b + 2) * (a + b + 3);
-      return  n / d;
+      using promoted_real_type = policies::evaluation_t<RealType, Policy>;
+      const auto a = static_cast<promoted_real_type>(dist.alpha());
+      const auto b = static_cast<promoted_real_type>(dist.beta());
+      const promoted_real_type a_2 = a * a;
+      const promoted_real_type n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2));
+      const promoted_real_type d = a * b * (a + b + 2) * (a + b + 3);
+      return static_cast<promoted_real_type>(n / d);
     } // kurtosis_excess
 
     template <class RealType, class Policy>
@@ -398,7 +409,7 @@ namespace boost
       {
         if (a == 1)
         {
-          return static_cast<RealType>(1 / beta(a, b));
+          return static_cast<RealType>(1 / beta(a, b, Policy()));
         }
         else if (a < 1)
         {
@@ -413,7 +424,7 @@ namespace boost
       {
         if (b == 1)
         {
-          return static_cast<RealType>(1 / beta(a, b));
+          return static_cast<RealType>(1 / beta(a, b, Policy()));
         }
         else if (b < 1)
         {

include/boost/math/policies/policy.hpp

@@ -766,6 +766,9 @@ struct evaluation<double, Policy>
    using type = typename boost::math::conditional<Policy::promote_double_type::value, long double, double>::type;
 };
 
+template <class Real, class Policy>
+using evaluation_t = typename evaluation<Real, Policy>::type;
+
 template <class Real, class Policy>
 struct precision
 {