diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 29feac48fe..4fe3c7d751 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -164,7 +164,104 @@ namespace boost return check_dist(function, x_min, x_max, result, pol) && check_prob(function, p, result, pol); } // bool check_dist_and_prob + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_mean(const RealType x_min, const RealType x_max) + { + return (x_min + x_max) / 2; + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_variance(const RealType x_min, const RealType x_max) + { + return (x_max - x_min) * (x_max - x_min) / 8; + } + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_pdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + using boost::math::constants::pi; + return 1 / (pi() * sqrt((x - x_min) * (x_max - x))); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_cdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (x == x_min) + { + return 0; + } + else if (x == x_max) + { + return 1; + } + using boost::math::constants::pi; + return 2 * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_ccdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases + if (x == x_min) + { + return 1; + } + else if (x == x_max) + { + return 0; + } + using boost::math::constants::pi; + // Naive version x = 1 - x; + // result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); + // is less accurate, so use acos instead of asin for complement. + return 2 * acos(sqrt((x - x_min) / (x_max - x_min))) / pi(); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_quantile( const RealType p, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (p == 0) + { + return 0; + } + if (p == 1) + { + return 1; + } + using boost::math::constants::half_pi; + RealType sin2hpip = sin(half_pi() * p); + RealType sin2hpip2 = sin2hpip * sin2hpip; + return -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; + } + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_cquantile(const RealType q, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (q == 1) + { + return 0; + } + if (q == 0) + { + return 1; + } + // Naive RealType p = 1 - q; result = sin(half_pi() * p); loses accuracy, so use a cos alternative instead. + //result = cos(half_pi() * q); // for arcsine(0,1) + //result = result * result; + // For generalized arcsine: + using boost::math::constants::half_pi; + RealType cos2hpip = cos(half_pi() * q); + RealType cos2hpip2 = cos2hpip * cos2hpip; + return -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; + } } // namespace arcsine_detail template > @@ -230,7 +327,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::mean(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -239,7 +336,9 @@ namespace boost { return result; } - return (x_min + x_max) / 2; + + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } // mean template @@ -248,7 +347,7 @@ namespace boost RealType result; RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::variance(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -257,7 +356,9 @@ namespace boost { return result; } - return (x_max - x_min) * (x_max - x_min) / 8; + + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_variance(static_cast(x_min), static_cast(x_max))); } // variance template @@ -277,7 +378,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); RealType result; - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::median(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -286,7 +387,9 @@ namespace boost { return result; } - return (x_min + x_max) / 2; + // The median is the same as the mean + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } template @@ -296,7 +399,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::skewness(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -315,7 +418,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::kurtosis_excess(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -324,8 +427,7 @@ namespace boost { return result; } - result = -3; - return result / 2; + return static_cast(-1.5f); } // kurtosis_excess template @@ -335,7 +437,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::kurtosis(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -344,75 +446,61 @@ namespace boost { return result; } - - return 3 + kurtosis_excess(dist); + return static_cast(1.5f); } // kurtosis template BOOST_MATH_GPU_ENABLED inline RealType pdf(const arcsine_distribution& dist, const RealType& xx) { // Probability Density/Mass Function arcsine. BOOST_FPU_EXCEPTION_GUARD - BOOST_MATH_STD_USING // For ADL of std functions. - - constexpr auto function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; - RealType lo = dist.x_min(); - RealType hi = dist.x_max(); + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); RealType x = xx; // Argument checks: - RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - lo, hi, x, + RealType result; + if (!arcsine_detail::check_dist_and_x( + "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - using boost::math::constants::pi; - result = static_cast(1) / (pi() * sqrt((x - lo) * (hi - x))); - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_pdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // pdf template BOOST_MATH_GPU_ENABLED inline RealType cdf(const arcsine_distribution& dist, const RealType& x) { // Cumulative Distribution Function arcsine. - BOOST_MATH_STD_USING // For ADL of std functions. - - constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; - RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); // Argument checks: RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - x_min, x_max, x, + if (!arcsine_detail::check_dist_and_x( + "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - // Special cases: - if (x == x_min) - { - return 0; - } - else if (x == x_max) - { - return 1; - } - using boost::math::constants::pi; - result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_cdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // arcsine cdf template BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type, RealType>& c) { // Complemented Cumulative Distribution Function arcsine. - BOOST_MATH_STD_USING // For ADL of std functions. - constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; - RealType x = c.param; arcsine_distribution const& dist = c.dist; RealType x_min = dist.x_min(); @@ -420,27 +508,19 @@ namespace boost // Argument checks: RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - x_min, x_max, x, + if (!arcsine_detail::check_dist_and_x( + "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - if (x == x_min) - { - return 0; - } - else if (x == x_max) - { - return 1; - } - using boost::math::constants::pi; - // Naive version x = 1 - x; - // result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); - // is less accurate, so use acos instead of asin for complement. - result = static_cast(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi(); - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_ccdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // arcsine ccdf template @@ -454,37 +534,22 @@ namespace boost // and return a value such that the probability that a random variable x // will be less than or equal to that value // is whatever probability you supplied as an argument. - BOOST_MATH_STD_USING // For ADL of std functions. - - using boost::math::constants::half_pi; - - constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; - RealType result = 0; // of argument checks: RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist_and_prob( - function, - x_min, x_max, p, + if (!arcsine_detail::check_dist_and_prob( + "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + p, &result, Policy())) { return result; } - // Special cases: - if (p == 0) - { - return 0; - } - if (p == 1) - { - return 1; - } - - RealType sin2hpip = sin(half_pi() * p); - RealType sin2hpip2 = sin2hpip * sin2hpip; - result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; - - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_quantile(static_cast(p), + static_cast(x_min), + static_cast(x_max))); } // quantile template @@ -493,19 +558,14 @@ namespace boost // Complement Quantile or Percent Point arcsine function. // Return the number of expected x for a given // complement of the probability q. - BOOST_MATH_STD_USING // For ADL of std functions. - - using boost::math::constants::half_pi; - constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; - // Error checks: RealType q = c.param; const arcsine_distribution& dist = c.dist; RealType result = 0; RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist_and_prob( - function, + if (!arcsine_detail::check_dist_and_prob( + "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)", x_min, x_max, q, @@ -513,26 +573,11 @@ namespace boost { return result; } - // Special cases: - if (q == 1) - { - return 0; - } - if (q == 0) - { - return 1; - } - // Naive RealType p = 1 - q; result = sin(half_pi() * p); loses accuracy, so use a cos alternative instead. - //result = cos(half_pi() * q); // for arcsine(0,1) - //result = result * result; - // For generalized arcsine: - RealType cos2hpip = cos(half_pi() * q); - RealType cos2hpip2 = cos2hpip * cos2hpip; - result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; - - return result; + typedef typename policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_cquantile(static_cast(q), + static_cast(x_min), + static_cast(x_max))); } // Quantile Complement - } // namespace math } // namespace boost diff --git a/include/boost/math/policies/policy.hpp b/include/boost/math/policies/policy.hpp index ccc33bb9f8..6d80b2b292 100644 --- a/include/boost/math/policies/policy.hpp +++ b/include/boost/math/policies/policy.hpp @@ -766,6 +766,9 @@ struct evaluation using type = typename boost::math::conditional::type; }; +template +using evaluation_t = typename evaluation::type; + template struct precision { diff --git a/test/test_arcsine.cpp b/test/test_arcsine.cpp index eed8dc9481..cb5646e627 100644 --- a/test/test_arcsine.cpp +++ b/test/test_arcsine.cpp @@ -437,11 +437,19 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.05), static_cast(0.85643370687129372924905811522494428117838480010259L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.5), static_cast(0.5L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.95), static_cast(0.14356629312870627075094188477505571882161519989741L), 8 * tolerance); // Not much less accurate. + BOOST_CHECK_EQUAL(cdf(as_m2m1, as_m2m1.x_min()), 0); + BOOST_CHECK_EQUAL(cdf(as_m2m1, as_m2m1.x_max()), 1); + BOOST_CHECK_EQUAL(cdf(complement(as_m2m1, as_m2m1.x_min())), 1); + BOOST_CHECK_EQUAL(cdf(complement(as_m2m1, as_m2m1.x_max())), 0); // Quantile BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.85643370687129372924905811522494428117838480010259L)), -static_cast(1.05L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.5L)), -static_cast(1.5L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.14356629312870627075094188477505571882161519989741L)), -static_cast(1.95L), 4 * tolerance); // + BOOST_CHECK_EQUAL(quantile(as_m2m1, 1), 1); + BOOST_CHECK_EQUAL(quantile(as_m2m1, 0), 0); + BOOST_CHECK_EQUAL(quantile(complement(as_m2m1, 1)), 0); + BOOST_CHECK_EQUAL(quantile(complement(as_m2m1, 0)), 1); BOOST_CHECK_CLOSE_FRACTION(quantile(complement(as_m2m1, static_cast(0.14356629312870627075094188477505571882161519989741L))), -static_cast(1.05L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.5L)), -static_cast(1.5L), 2 * tolerance); //