diff --git a/include/boost/math/ccmath/abs.hpp b/include/boost/math/ccmath/abs.hpp
index fa798eba12..d8077fc3c5 100644
--- a/include/boost/math/ccmath/abs.hpp
+++ b/include/boost/math/ccmath/abs.hpp
@@ -12,6 +12,7 @@
 #include <type_traits>
 #include <limits>
 #include <boost/math/tools/is_constant_evaluated.hpp>
+#include <boost/math/concepts/concepts.hpp>
 #include <boost/math/tools/assert.hpp>
 #include <boost/math/ccmath/isnan.hpp>
 #include <boost/math/ccmath/isinf.hpp>
@@ -28,8 +29,8 @@ namespace boost::math::ccmath {
 
 namespace detail {
 
-template <typename T> 
-constexpr T abs_impl(T x) noexcept
+template <BOOST_MATH_ARBITRARY_SIGNED_ARITHMETIC T> 
+inline constexpr T abs_impl(T x) noexcept
 {
     if (boost::math::ccmath::isnan(x))
     {
@@ -51,7 +52,8 @@ constexpr T abs_impl(T x) noexcept
 } // Namespace detail
 
 template <typename T, std::enable_if_t<!std::is_unsigned_v<T>, bool> = true>
-constexpr T abs(T x) noexcept
+    BOOST_MATH_REQUIRES(BOOST_MATH_ARBITRARY_SIGNED_ARITHMETIC, T)
+inline constexpr T abs(T x) noexcept
 {
     if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
     {
@@ -67,7 +69,8 @@ constexpr T abs(T x) noexcept
 // If abs() is called with an argument of type X for which is_unsigned_v<X> is true and if X
 // cannot be converted to int by integral promotion (7.3.7), the program is ill-formed.
 template <typename T, std::enable_if_t<std::is_unsigned_v<T>, bool> = true>
-constexpr T abs(T x) noexcept
+    BOOST_MATH_REQUIRES(BOOST_MATH_ARBITRARY_UNSIGNED_ARITHMETIC, T)
+inline constexpr T abs(T x) noexcept
 {
     if constexpr (std::is_convertible_v<T, int>)
     {
diff --git a/include/boost/math/concepts/concepts.hpp b/include/boost/math/concepts/concepts.hpp
new file mode 100644
index 0000000000..e0ee3646ea
--- /dev/null
+++ b/include/boost/math/concepts/concepts.hpp
@@ -0,0 +1,378 @@
+//  (C) Copyright Matt Borland 2022.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_CONCEPTS_CONCEPTS_HPP
+#define BOOST_MATH_CONCEPTS_CONCEPTS_HPP
+
+#if (__cplusplus >= 202002L || _MSVC_LANG >= 202002L) && !defined(BOOST_MATH_DISABLE_CONCEPTS)
+#if __has_include(<concepts>)
+
+#include <utility>
+#include <algorithm>
+#include <concepts>
+#include <functional>
+#include <type_traits>
+#include <limits>
+#include <iterator>
+#include <complex>
+#include <boost/math/tools/config.hpp>
+#include <boost/math/policies/policy.hpp>
+
+namespace boost::math::concepts {
+
+namespace detail {
+
+template <typename X, typename Y, typename Op>
+struct op_valid_impl
+{
+    template <typename U, typename L, typename R>
+    static constexpr auto test(int) -> decltype(std::declval<U>()(std::declval<L>(), std::declval<R>()),
+                                                void(), std::true_type());
+
+    template <typename U, typename L, typename R>
+    static constexpr auto test(...) -> std::false_type;
+
+    using type = decltype(test<Op, X, Y>(0));
+};
+
+template <typename X, typename Y, typename Op> 
+using op_valid_t = typename op_valid_impl<X, Y, Op>::type;
+
+template <typename X, typename Y, typename Op>
+inline constexpr bool op_valid_v = op_valid_t<X, Y, Op>::value;
+
+// Detector for class member functions
+struct nonesuch 
+{
+    nonesuch(const nonesuch&) = delete;
+    ~nonesuch() = delete;
+    void operator=(const nonesuch&) = delete;
+};
+
+template <typename Default, typename AlwaysVoid, template<typename...> typename Op, typename... Args>
+struct detector 
+{
+    using value_t = std::false_type;
+    using type = Default;
+};
+ 
+template <typename Default, template<typename...> typename Op, typename... Args>
+struct detector<Default, std::void_t<Op<Args...>>, Op, Args...> 
+{
+    using value_t = std::true_type;
+    using type = Op<Args...>;
+};
+
+template <template<typename...> typename Op, typename... Args>
+using is_detected = typename detector<nonesuch, void, Op, Args...>::value_t;
+ 
+template <template<typename...> typename Op, typename... Args>
+using detected_t = typename detector<nonesuch, void, Op, Args...>::type;
+
+#define BOOST_MATH_HAS_MEMBER_FUNCTION(member)                                      \
+template <typename T>                                                               \
+using has_##member##_t = decltype(std::declval<T&>().member());                     \
+                                                                                    \
+template <typename T>                                                               \
+inline constexpr bool has_##member##_v = is_detected<has_##member##_t, T>::value;       
+
+BOOST_MATH_HAS_MEMBER_FUNCTION(begin)
+BOOST_MATH_HAS_MEMBER_FUNCTION(end)
+BOOST_MATH_HAS_MEMBER_FUNCTION(real)
+BOOST_MATH_HAS_MEMBER_FUNCTION(imag)
+
+} // Namespace detail
+
+template <typename T>
+concept integral = std::is_integral_v<T>;
+
+template <typename T>
+concept signed_integral = integral<T> && std::is_signed_v<T>;
+
+template <typename T>
+concept unsigned_integral = integral<T> && std::is_unsigned_v<T>;
+
+template <typename T>
+concept real = std::is_floating_point_v<T>;
+
+template <typename T>
+concept complex = std::is_same_v<T, std::complex<float>>
+                  || std::is_same_v<T, std::complex<double>>
+                  #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+                  || std::is_same_v<T, std::complex<long double>>
+                  #endif
+                  ;
+
+template <typename T>
+concept real_or_complex = real<T> || complex<T>;
+
+template <typename T>
+concept arithmetic = std::is_arithmetic_v<T>;
+
+template <typename T>
+concept numerical = arithmetic<T> || real_or_complex<T>;
+
+template <typename T>
+concept signed_arithmetic = arithmetic<T> && std::is_signed_v<T>;
+
+template <typename T>
+concept unsigned_arithmetic = arithmetic<T> && std::is_unsigned_v<T>;
+
+template <typename T>
+concept arbitrary_unsigned_arithmetic_type = unsigned_arithmetic<T> ||
+                                             (detail::op_valid_v<T, T, std::equal_to<>> &&
+                                              detail::op_valid_v<T, T, std::not_equal_to<>> &&
+                                              detail::op_valid_v<T, T, std::greater<>> &&
+                                              detail::op_valid_v<T, T, std::less<>> &&
+                                              detail::op_valid_v<T, T, std::greater_equal<>> &&
+                                              detail::op_valid_v<T, T, std::less_equal<>> &&
+                                              detail::op_valid_v<T, T, std::plus<>> &&
+                                              detail::op_valid_v<T, T, std::minus<>> &&
+                                              detail::op_valid_v<T, T, std::multiplies<>> &&
+                                              detail::op_valid_v<T, T, std::divides<>>);
+
+template <typename T>
+concept arbitrary_signed_arithmetic_type = signed_arithmetic<T> ||
+                                           (arbitrary_unsigned_arithmetic_type<T> &&
+                                            (detail::op_valid_v<T, T, std::negate<>> ||
+                                             std::numeric_limits<T>::is_signed));
+
+template <typename T>
+concept arbitrary_arithmetic_type = arbitrary_unsigned_arithmetic_type<T> ||
+                                    arbitrary_signed_arithmetic_type<T>;
+
+template <typename T>
+concept arbitrary_unsigned_integer_type = arbitrary_unsigned_arithmetic_type<T> &&
+                                           std::numeric_limits<T>::is_integer;
+
+template <typename T>
+concept arbitrary_signed_integer_type = arbitrary_signed_arithmetic_type<T> &&
+                                         std::numeric_limits<T>::is_integer;
+
+template <typename T>
+concept arbitrary_integer_type = arbitrary_unsigned_integer_type<T> ||
+                                  arbitrary_signed_integer_type<T>;
+
+template <typename T>
+concept arbitrary_real_type = arbitrary_arithmetic_type<T> &&
+                               !std::numeric_limits<T>::is_integer;
+
+template <typename T>
+concept arbitrary_complex_type = complex<T> ||
+                                 (detail::has_real_v<T> &&
+                                  detail::has_imag_v<T>);
+
+template <typename T>
+concept arbitrary_real_or_complex_type = arbitrary_real_type<T> ||
+                                         arbitrary_complex_type<T>;
+
+template <typename T>
+concept arbitrary_numerical_type = arbitrary_real_or_complex_type<T> ||
+                                   arbitrary_arithmetic_type<T>;
+
+template <typename T>
+concept policy = boost::math::policies::is_policy<T>::value;
+
+template <typename Derived, typename Base>
+concept derived_from = std::is_base_of_v<Base, Derived> &&
+                       std::is_convertible_v<const volatile Derived*, const volatile Base*>;
+
+template <typename T>
+concept forward_iterator = derived_from<typename std::iterator_traits<T>::iterator_category, std::forward_iterator_tag>;
+
+template <typename T>
+concept bidirectional_iterator = derived_from<typename std::iterator_traits<T>::iterator_category, std::bidirectional_iterator_tag>;
+
+template <typename T>
+concept random_access_iterator = derived_from<typename std::iterator_traits<T>::iterator_category, std::random_access_iterator_tag>;
+
+template <typename I, typename T>
+concept output_iterator = derived_from<typename std::iterator_traits<I>::iterator_category, std::input_iterator_tag> &&
+                          derived_from<typename std::iterator_traits<T>::iterator_category, std::output_iterator_tag>;
+
+template <typename T>
+concept is_container = detail::has_begin_v<T> &&
+                       detail::has_end_v<T>;
+
+template <typename T>
+concept random_access_container = is_container<T> &&
+                                  boost::math::concepts::random_access_iterator<typename T::iterator>;
+
+} // boost::math::concepts
+
+#define BOOST_MATH_INTEGRAL boost::math::concepts::integral
+#define BOOST_MATH_SIGNED_INTEGRAL boost::math::concepts::signed_integral
+#define BOOST_MATH_UNSIGNED_INTEGRAL boost::math::concepts::unsigned_integral
+#define BOOST_MATH_REAL boost::math::concepts::real
+#define BOOST_MATH_COMPLEX boost::math::concepts::complex
+#define BOOST_MATH_REAL_OR_COMPLEX boost::math::concepts::real_or_complex
+#define BOOST_MATH_ARITHMETIC boost::math::concepts::arithmetic
+#define BOOST_MATH_NUMERICAL boost::math::concepts::numerical
+#define BOOST_MATH_SIGNED_ARITHMETIC boost::math::concepts::signed_arithmetic
+#define BOOST_MATH_UNSIGNED_ARITHMETIC boost::math::concepts::unsigned_arithmetic
+#define BOOST_MATH_ARBITRARY_UNSIGNED_ARITHMETIC boost::math::concepts::arbitrary_unsigned_arithmetic_type
+#define BOOST_MATH_ARBITRARY_SIGNED_ARITHMETIC boost::math::concepts::arbitrary_signed_arithmetic_type
+#define BOOST_MATH_ARBITRARY_ARITHMETIC boost::math::concepts::arbitrary_arithmetic_type
+#define BOOST_MATH_ARBITRARY_UNSIGNED_INTEGER boost::math::concepts::arbitrary_unsigned_integer_type
+#define BOOST_MATH_ARBITRARY_SIGNED_INTEGER boost::math::concepts::arbitrary_signed_integer_type
+#define BOOST_MATH_ARBITRARY_INTEGER boost::math::concepts::arbitrary_integer_type
+#define BOOST_MATH_ARBITRARY_REAL boost::math::concepts::arbitrary_real_type
+#define BOOST_MATH_ARBITRARY_COMPLEX boost::math::concepts::arbitrary_complex_type
+#define BOOST_MATH_ARBITRARY_REAL_OR_COMPLEX boost::math::concepts::arbitrary_real_or_complex_type
+#define BOOST_MATH_ARBITRARY_NUMERICAL boost::math::concepts::arbitrary_numerical_type
+
+#define BOOST_MATH_POLICY boost::math::concepts::policy
+
+#define BOOST_MATH_CONTAINER boost::math::concepts::is_container
+#define BOOST_MATH_RANDOM_ACCESS_CONTAINER boost::math::concepts::random_access_container
+
+#define BOOST_MATH_FORWARD_ITER boost::math::concepts::forward_iterator
+#define BOOST_MATH_BIDIRECTIONAL_ITER boost::math::concepts::bidirectional_iterator
+#define BOOST_MATH_RANDOM_ACCESS_ITER boost::math::concepts::random_access_iterator
+#define BOOST_MATH_OUTPUT_ITER(I, T) boost::math::concepts::output_iterator<I, T>
+#define BOOST_MATH_REQUIRES_ITER(X) requires X
+
+#define BOOST_MATH_REQUIRES(X, T) requires X<T>
+
+#ifdef BOOST_MATH_EXEC_COMPATIBLE
+#include <execution>
+
+namespace boost::math::concepts {
+
+template <typename T>
+concept execution_policy = std::is_execution_policy_v<std::remove_cvref_t<T>>;
+
+} // Namespace boost::math::concepts
+
+#define BOOST_MATH_EXECUTION_POLICY boost::math::concepts::execution_policy
+
+#endif // Has <execution>
+
+#endif // Has <concepts>
+#endif // C++20
+
+// If concepts are unavailable replace them with typename for compatibility
+
+#ifndef BOOST_MATH_INTEGRAL
+#  define BOOST_MATH_INTEGRAL typename
+#endif
+
+#ifndef BOOST_MATH_SIGNED_INTEGRAL
+#  define BOOST_MATH_SIGNED_INTEGRAL typename
+#endif
+
+#ifndef BOOST_MATH_UNSIGNED_INTEGRAL
+#  define BOOST_MATH_UNSIGNED_INTEGRAL typename
+#endif
+
+#ifndef BOOST_MATH_REAL
+#  define BOOST_MATH_REAL typename
+#endif
+
+#ifndef BOOST_MATH_COMPLEX
+#  define BOOST_MATH_COMPLEX typename
+#endif
+
+#ifndef BOOST_MATH_REAL_OR_COMPLEX
+#  define BOOST_MATH_REAL_OR_COMPLEX typename
+#endif
+
+#ifndef BOOST_MATH_ARITHMETIC
+#  define BOOST_MATH_ARITHMETIC typename
+#endif
+
+#ifndef BOOST_MATH_NUMERICAL
+#  define BOOST_MATH_NUMERICAL typename
+#endif
+
+#ifndef BOOST_MATH_SIGNED_ARITHMETIC
+#  define BOOST_MATH_SIGNED_ARITHMETIC typename
+#endif
+
+#ifndef BOOST_MATH_UNSIGNED_ARITHMETIC
+#  define BOOST_MATH_UNSIGNED_ARITHMETIC typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_UNSIGNED_ARITHMETIC
+#  define BOOST_MATH_ARBITRARY_UNSIGNED_ARITHMETIC typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_SIGNED_ARITHMETIC
+#  define BOOST_MATH_ARBITRARY_SIGNED_ARITHMETIC typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_ARITHMETIC
+#  define BOOST_MATH_ARBITRARY_ARITHMETIC typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_UNSIGNED_INTEGER
+#  define BOOST_MATH_ARBITRARY_UNSIGNED_INTEGER typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_SIGNED_INTEGER
+#  define BOOST_MATH_ARBITRARY_SIGNED_INTEGER typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_INTEGER
+#  define BOOST_MATH_ARBITRARY_INTEGER typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_REAL
+#  define BOOST_MATH_ARBITRARY_REAL typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_COMPLEX
+#  define BOOST_MATH_ARBITRARY_COMPLEX typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_REAL_OR_COMPLEX
+#  define BOOST_MATH_ARBITRARY_REAL_OR_COMPLEX typename
+#endif
+
+#ifndef BOOST_MATH_ARBITRARY_NUMERICAL
+#  define BOOST_MATH_ARBITRARY_NUMERICAL typename
+#endif
+
+#ifndef BOOST_MATH_POLICY
+#  define BOOST_MATH_POLICY typename
+#endif
+
+#ifndef BOOST_MATH_FORWARD_ITER
+#  define BOOST_MATH_FORWARD_ITER typename
+#endif
+
+#ifndef BOOST_MATH_BIDIRECTIONAL_ITER
+#  define BOOST_MATH_BIDIRECTIONAL_ITER typename
+#endif
+
+#ifndef BOOST_MATH_RANDOM_ACCESS_ITER
+#  define BOOST_MATH_RANDOM_ACCESS_ITER typename
+#endif
+
+#ifndef BOOST_MATH_OUTPUT_ITER
+#  define BOOST_MATH_OUTPUT_ITER(I, T)
+#endif
+
+#ifndef BOOST_MATH_REQUIRES_ITER
+#  define BOOST_MATH_REQUIRES_ITER(X)
+#endif
+
+#ifndef BOOST_MATH_CONTAINER
+#  define BOOST_MATH_CONTAINER typename
+#endif
+
+#ifndef BOOST_MATH_RANDOM_ACCESS_CONTAINER
+#  define BOOST_MATH_RANDOM_ACCESS_CONTAINER typename
+#endif
+
+#ifndef BOOST_MATH_EXECUTION_POLICY
+#  define BOOST_MATH_EXECUTION_POLICY typename
+#endif
+
+#ifndef BOOST_MATH_REQUIRES
+#  define BOOST_MATH_REQUIRES(X, T)
+#endif
+
+#endif // BOOST_MATH_CONCEPTS_CONCEPTS_HPP
diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index c950670383..b56358c2c5 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/include/boost/math/special_functions/beta.hpp
@@ -21,6 +21,8 @@
 #include <boost/math/special_functions/trunc.hpp>
 #include <boost/math/tools/roots.hpp>
 #include <boost/math/tools/assert.hpp>
+#include <boost/math/concepts/concepts.hpp>
+#include <limits>
 #include <cmath>
 
 namespace boost{ namespace math{
@@ -30,7 +32,7 @@ namespace detail{
 //
 // Implementation of Beta(a,b) using the Lanczos approximation:
 //
-template <class T, class Lanczos, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, class Lanczos, BOOST_MATH_POLICY Policy>
 T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
 {
    BOOST_MATH_STD_USING  // for ADL of std names
@@ -120,7 +122,7 @@ T beta_imp(T a, T b, const Lanczos&, const Policy& pol)
 // Generic implementation of Beta(a,b) without Lanczos approximation support
 // (Caution this is slow!!!):
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T beta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol)
 {
    BOOST_MATH_STD_USING
@@ -203,7 +205,7 @@ T beta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol)
 // powers are *hard* though, and using logarithms just leads to
 // horrendous cancellation errors.
 //
-template <class T, class Lanczos, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, class Lanczos, BOOST_MATH_POLICY Policy>
 T ibeta_power_terms(T a,
                         T b,
                         T x,
@@ -440,7 +442,7 @@ T ibeta_power_terms(T a,
 //
 // This version is generic, slow, and does not use the Lanczos approximation.
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T ibeta_power_terms(T a,
                         T b,
                         T x,
@@ -532,7 +534,7 @@ T ibeta_power_terms(T a,
 //
 // Series approximation to the incomplete beta:
 //
-template <class T>
+template <BOOST_MATH_ARBITRARY_REAL T>
 struct ibeta_series_t
 {
    typedef T result_type;
@@ -551,7 +553,7 @@ struct ibeta_series_t
    int n;
 };
 
-template <class T, class Lanczos, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, class Lanczos, BOOST_MATH_POLICY Policy>
 T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol)
 {
    BOOST_MATH_STD_USING
@@ -620,7 +622,7 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
 //
 // Incomplete Beta series again, this time without Lanczos support:
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos& l, bool normalised, T* p_derivative, T y, const Policy& pol)
 {
    BOOST_MATH_STD_USING
@@ -682,7 +684,7 @@ T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczo
 //
 // Continued fraction for the incomplete beta:
 //
-template <class T>
+template <BOOST_MATH_ARBITRARY_REAL T>
 struct ibeta_fraction2_t
 {
    typedef std::pair<T, T> result_type;
@@ -711,7 +713,7 @@ struct ibeta_fraction2_t
 //
 // Evaluate the incomplete beta via the continued fraction representation:
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative)
 {
    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
@@ -734,7 +736,7 @@ inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised,
 //
 // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x):
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative)
 {
    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
@@ -768,7 +770,7 @@ T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T*
 // beta(a,b,x) = prefix + delta * beta(a+k,b,x)
 // it is currently only called for small k.
 //
-template <class T>
+template <BOOST_MATH_ARBITRARY_REAL T>
 inline T rising_factorial_ratio(T a, T b, int k)
 {
    // calculate:
@@ -795,7 +797,7 @@ inline T rising_factorial_ratio(T a, T b, int k)
 // Note that the table size should never exceed the size of our
 // tables of factorials.
 //
-template <class T>
+template <BOOST_MATH_ARBITRARY_REAL T>
 struct Pn_size
 {
    // This is likely to be enough for ~35-50 digit accuracy
@@ -825,7 +827,7 @@ struct Pn_size<long double>
    static_assert(::boost::math::max_factorial<long double>::value >= 100, "Type does not provide for ~35-50 digits of accuracy");
 };
 
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised)
 {
    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
@@ -937,7 +939,7 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& po
 // For integer arguments we can relate the incomplete beta to the
 // complement of the binomial distribution cdf and use this finite sum.
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T binomial_ccdf(T n, T k, T x, T y, const Policy& pol)
 {
    BOOST_MATH_STD_USING // ADL of std names
@@ -998,7 +1000,7 @@ T binomial_ccdf(T n, T k, T x, T y, const Policy& pol)
 // input range and select the right implementation method for
 // each domain:
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative)
 {
    static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
@@ -1398,13 +1400,13 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de
    return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
 } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)
 
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)
 {
    return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(nullptr));
 }
 
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_REAL T, BOOST_MATH_POLICY Policy>
 T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
 {
    static const char* function = "ibeta_derivative<%1%>(%1%,%1%,%1%)";
diff --git a/include/boost/math/special_functions/log1p.hpp b/include/boost/math/special_functions/log1p.hpp
index 8121a573cd..877ec214d9 100644
--- a/include/boost/math/special_functions/log1p.hpp
+++ b/include/boost/math/special_functions/log1p.hpp
@@ -23,6 +23,7 @@
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/tools/assert.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/concepts/concepts.hpp>
 
 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
 //
@@ -293,7 +294,7 @@ const typename log1p_initializer<T, Policy, tag>::init log1p_initializer<T, Poli
 
 } // namespace detail
 
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_ARITHMETIC T, BOOST_MATH_POLICY Policy>
 inline typename tools::promote_args<T>::type log1p(T x, const Policy&)
 {
    typedef typename tools::promote_args<T>::type result_type;
@@ -327,7 +328,7 @@ inline typename tools::promote_args<T>::type log1p(T x, const Policy&)
 
 #if defined(BOOST_HAS_LOG1P) && !(defined(__osf__) && defined(__DECCXX_VER))
 #  ifdef BOOST_MATH_USE_C99
-template <class Policy>
+template <BOOST_MATH_POLICY Policy>
 inline float log1p(float x, const Policy& pol)
 {
    if(x < -1)
@@ -339,7 +340,7 @@ inline float log1p(float x, const Policy& pol)
    return ::log1pf(x);
 }
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-template <class Policy>
+template <BOOST_MATH_POLICY Policy>
 inline long double log1p(long double x, const Policy& pol)
 {
    if(x < -1)
@@ -352,7 +353,7 @@ inline long double log1p(long double x, const Policy& pol)
 }
 #endif
 #else
-template <class Policy>
+template <BOOST_MATH_POLICY Policy>
 inline float log1p(float x, const Policy& pol)
 {
    if(x < -1)
@@ -364,7 +365,7 @@ inline float log1p(float x, const Policy& pol)
    return ::log1p(x);
 }
 #endif
-template <class Policy>
+template <BOOST_MATH_POLICY Policy>
 inline double log1p(double x, const Policy& pol)
 {
    if(x < -1)
@@ -381,7 +382,7 @@ inline double log1p(double x, const Policy& pol)
 // that your compilers optimizer won't mess this code up!!
 // Currently tested with VC8 and Intel 9.1.
 //
-template <class Policy>
+template <BOOST_MATH_POLICY Policy>
 inline double log1p(double x, const Policy& pol)
 {
    if(x < -1)
@@ -396,7 +397,7 @@ inline double log1p(double x, const Policy& pol)
    else
       return ::log(u)*(x/(u-1.0));
 }
-template <class Policy>
+template <BOOST_MATH_POLICY Policy>
 inline float log1p(float x, const Policy& pol)
 {
    return static_cast<float>(boost::math::log1p(static_cast<double>(x), pol));
@@ -406,7 +407,7 @@ inline float log1p(float x, const Policy& pol)
 // For some reason this fails to compile under WinCE...
 // Needs more investigation.
 //
-template <class Policy>
+template <BOOST_MATH_POLICY Policy>
 inline long double log1p(long double x, const Policy& pol)
 {
    if(x < -1)
@@ -424,7 +425,7 @@ inline long double log1p(long double x, const Policy& pol)
 #endif
 #endif
 
-template <class T>
+template <BOOST_MATH_ARBITRARY_ARITHMETIC T>
 inline typename tools::promote_args<T>::type log1p(T x)
 {
    return boost::math::log1p(x, policies::policy<>());
@@ -432,7 +433,7 @@ inline typename tools::promote_args<T>::type log1p(T x)
 //
 // Compute log(1+x)-x:
 //
-template <class T, class Policy>
+template <BOOST_MATH_ARBITRARY_ARITHMETIC T, BOOST_MATH_POLICY Policy>
 inline typename tools::promote_args<T>::type
    log1pmx(T x, const Policy& pol)
 {
@@ -464,7 +465,7 @@ inline typename tools::promote_args<T>::type
    return result;
 }
 
-template <class T>
+template <BOOST_MATH_ARBITRARY_ARITHMETIC T>
 inline typename tools::promote_args<T>::type log1pmx(T x)
 {
    return log1pmx(x, policies::policy<>());
diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp
index 6f9d739e1e..0ba43b5ae5 100644
--- a/include/boost/math/special_functions/math_fwd.hpp
+++ b/include/boost/math/special_functions/math_fwd.hpp
@@ -29,6 +29,7 @@
 #include <boost/math/special_functions/detail/round_fwd.hpp>
 #include <boost/math/tools/promotion.hpp> // for argument promotion.
 #include <boost/math/policies/policy.hpp>
+#include <boost/math/concepts/concepts.hpp>
 
 #define BOOST_NO_MACRO_EXPAND /**/
 
@@ -567,17 +568,17 @@ namespace boost
    typename tools::promote_args<RT>::type cbrt(RT z, const Policy&);
 
    // log1p is log(x + 1)
-   template <class T>
+   template <BOOST_MATH_ARBITRARY_ARITHMETIC T>
    typename tools::promote_args<T>::type log1p(T);
 
-   template <class T, class Policy>
+   template <BOOST_MATH_ARBITRARY_ARITHMETIC T, BOOST_MATH_POLICY Policy>
    typename tools::promote_args<T>::type log1p(T, const Policy&);
 
    // log1pmx is log(x + 1) - x
-   template <class T>
+   template <BOOST_MATH_ARBITRARY_ARITHMETIC T>
    typename tools::promote_args<T>::type log1pmx(T);
 
-   template <class T, class Policy>
+   template <BOOST_MATH_ARBITRARY_ARITHMETIC T, BOOST_MATH_POLICY Policy>
    typename tools::promote_args<T>::type log1pmx(T, const Policy&);
 
    // Exp (x) minus 1 functions.
diff --git a/include/boost/math/statistics/univariate_statistics.hpp b/include/boost/math/statistics/univariate_statistics.hpp
index 082517d650..35a0b25ac6 100644
--- a/include/boost/math/statistics/univariate_statistics.hpp
+++ b/include/boost/math/statistics/univariate_statistics.hpp
@@ -10,6 +10,7 @@
 #include <boost/math/statistics/detail/single_pass.hpp>
 #include <boost/math/tools/config.hpp>
 #include <boost/math/tools/assert.hpp>
+#include <boost/math/concepts/concepts.hpp>
 #include <algorithm>
 #include <iterator>
 #include <tuple>
@@ -25,7 +26,7 @@
 
 namespace boost::math::statistics {
 
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto mean(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
@@ -55,25 +56,25 @@ inline auto mean(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator
     }
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto mean(ExecutionPolicy&& exec, Container const & v)
 {
     return mean(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto mean(ForwardIterator first, ForwardIterator last)
 {
     return mean(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto mean(Container const & v)
 {
     return mean(std::execution::seq, std::cbegin(v), std::cend(v));
 }
 
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
@@ -104,25 +105,25 @@ inline auto variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardItera
     }
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto variance(ExecutionPolicy&& exec, Container const & v)
 {
     return variance(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto variance(ForwardIterator first, ForwardIterator last)
 {
     return variance(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto variance(Container const & v)
 {
     return variance(std::execution::seq, std::cbegin(v), std::cend(v));
 }
 
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto sample_variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     const auto n = std::distance(first, last);
@@ -130,25 +131,25 @@ inline auto sample_variance(ExecutionPolicy&& exec, ForwardIterator first, Forwa
     return n*variance(exec, first, last)/(n-1);
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto sample_variance(ExecutionPolicy&& exec, Container const & v)
 {
     return sample_variance(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto sample_variance(ForwardIterator first, ForwardIterator last)
 {
     return sample_variance(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto sample_variance(Container const & v)
 {
     return sample_variance(std::execution::seq, std::cbegin(v), std::cend(v));
 }
 
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto mean_and_sample_variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
@@ -181,25 +182,25 @@ inline auto mean_and_sample_variance(ExecutionPolicy&& exec, ForwardIterator fir
     }
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto mean_and_sample_variance(ExecutionPolicy&& exec, Container const & v)
 {
     return mean_and_sample_variance(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto mean_and_sample_variance(ForwardIterator first, ForwardIterator last)
 {
     return mean_and_sample_variance(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto mean_and_sample_variance(Container const & v)
 {
     return mean_and_sample_variance(std::execution::seq, std::cbegin(v), std::cend(v));
 }
 
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto first_four_moments(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
@@ -236,26 +237,26 @@ inline auto first_four_moments(ExecutionPolicy&& exec, ForwardIterator first, Fo
     }
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto first_four_moments(ExecutionPolicy&& exec, Container const & v)
 {
     return first_four_moments(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto first_four_moments(ForwardIterator first, ForwardIterator last)
 {
     return first_four_moments(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto first_four_moments(Container const & v)
 {
     return first_four_moments(std::execution::seq, std::cbegin(v), std::cend(v));
 }
 
 // https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto skewness(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
@@ -298,19 +299,19 @@ inline auto skewness(ExecutionPolicy&& exec, ForwardIterator first, ForwardItera
     }
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto skewness(ExecutionPolicy&& exec, Container & v)
 {
     return skewness(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto skewness(ForwardIterator first, ForwardIterator last)
 {
     return skewness(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto skewness(Container const & v)
 {
     return skewness(std::execution::seq, std::cbegin(v), std::cend(v));
@@ -318,7 +319,7 @@ inline auto skewness(Container const & v)
 
 // Follows equation 1.6 of:
 // https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto kurtosis(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     const auto [M1, M2, M3, M4] = first_four_moments(exec, first, last);
@@ -329,50 +330,50 @@ inline auto kurtosis(ExecutionPolicy&& exec, ForwardIterator first, ForwardItera
     return M4/(M2*M2);
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto kurtosis(ExecutionPolicy&& exec, Container const & v)
 {
     return kurtosis(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto kurtosis(ForwardIterator first, ForwardIterator last)
 {
     return kurtosis(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto kurtosis(Container const & v)
 {
     return kurtosis(std::execution::seq, std::cbegin(v), std::cend(v));
 }
 
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto excess_kurtosis(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     return kurtosis(exec, first, last) - 3;
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto excess_kurtosis(ExecutionPolicy&& exec, Container const & v)
 {
     return excess_kurtosis(exec, std::cbegin(v), std::cend(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto excess_kurtosis(ForwardIterator first, ForwardIterator last)
 {
     return excess_kurtosis(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto excess_kurtosis(Container const & v)
 {
     return excess_kurtosis(std::execution::seq, std::cbegin(v), std::cend(v));
 }
 
 
-template<class ExecutionPolicy, class RandomAccessIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 auto median(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last)
 {
     const auto num_elems = std::distance(first, last);
@@ -393,25 +394,25 @@ auto median(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIter
 }
 
 
-template<class ExecutionPolicy, class RandomAccessContainer>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto median(ExecutionPolicy&& exec, RandomAccessContainer & v)
 {
     return median(exec, std::begin(v), std::end(v));
 }
 
-template<class RandomAccessIterator>
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 inline auto median(RandomAccessIterator first, RandomAccessIterator last)
 {
     return median(std::execution::seq, first, last);
 }
 
-template<class RandomAccessContainer>
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto median(RandomAccessContainer & v)
 {
     return median(std::execution::seq, std::begin(v), std::end(v));
 }
 
-template<class ExecutionPolicy, class RandomAccessIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last)
 {
     using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
@@ -444,50 +445,50 @@ inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first,
     }
 }
 
-template<class ExecutionPolicy, class RandomAccessContainer>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessContainer & v)
 {
     return gini_coefficient(exec, std::begin(v), std::end(v));
 }
 
-template<class RandomAccessIterator>
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 inline auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 {
     return gini_coefficient(std::execution::seq, first, last);
 }
 
-template<class RandomAccessContainer>
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto gini_coefficient(RandomAccessContainer & v)
 {
     return gini_coefficient(std::execution::seq, std::begin(v), std::end(v));
 }
 
-template<class ExecutionPolicy, class RandomAccessIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 inline auto sample_gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last)
 {
     const auto n = std::distance(first, last);
     return n*gini_coefficient(exec, first, last)/(n-1);
 }
 
-template<class ExecutionPolicy, class RandomAccessContainer>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto sample_gini_coefficient(ExecutionPolicy&& exec, RandomAccessContainer & v)
 {
     return sample_gini_coefficient(exec, std::begin(v), std::end(v));
 }
 
-template<class RandomAccessIterator>
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 {
     return sample_gini_coefficient(std::execution::seq, first, last);
 }
 
-template<class RandomAccessContainer>
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto sample_gini_coefficient(RandomAccessContainer & v)
 {
     return sample_gini_coefficient(std::execution::seq, std::begin(v), std::end(v));
 }
 
-template<class ExecutionPolicy, class RandomAccessIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last,
     typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())
 {
@@ -516,28 +517,28 @@ auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessIterator firs
     }
 }
 
-template<class ExecutionPolicy, class RandomAccessContainer>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessContainer & v,
     typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())
 {
     return median_absolute_deviation(exec, std::begin(v), std::end(v), center);
 }
 
-template<class RandomAccessIterator>
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 inline auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last,
     typename RandomAccessIterator::value_type center=std::numeric_limits<typename RandomAccessIterator::value_type>::quiet_NaN())
 {
     return median_absolute_deviation(std::execution::seq, first, last, center);
 }
 
-template<class RandomAccessContainer>
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto median_absolute_deviation(RandomAccessContainer & v,
     typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())
 {
     return median_absolute_deviation(std::execution::seq, std::begin(v), std::end(v), center);
 }
 
-template<class ExecutionPolicy, class ForwardIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator>
 auto interquartile_range(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
@@ -578,25 +579,26 @@ auto interquartile_range(ExecutionPolicy&& exec, ForwardIterator first, ForwardI
     }
 }
 
-template<class ExecutionPolicy, class RandomAccessContainer>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto interquartile_range(ExecutionPolicy&& exec, RandomAccessContainer & v)
 {
     return interquartile_range(exec, std::begin(v), std::end(v));
 }
 
-template<class RandomAccessIterator>
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator>
 inline auto interquartile_range(RandomAccessIterator first, RandomAccessIterator last)
 {
     return interquartile_range(std::execution::seq, first, last);
 }
 
-template<class RandomAccessContainer>
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer>
 inline auto interquartile_range(RandomAccessContainer & v)
 {
     return interquartile_range(std::execution::seq, std::begin(v), std::end(v));
 }
 
-template<class ExecutionPolicy, class ForwardIterator, class OutputIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator, class OutputIterator>
+    BOOST_MATH_REQUIRES_ITER(BOOST_MATH_OUTPUT_ITER(ForwardIterator, OutputIterator))
 inline OutputIterator mode(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last, OutputIterator output)
 {
     if(!std::is_sorted(exec, first, last))
@@ -614,13 +616,15 @@ inline OutputIterator mode(ExecutionPolicy&& exec, ForwardIterator first, Forwar
     return detail::mode_impl(first, last, output);
 }
 
-template<class ExecutionPolicy, class Container, class OutputIterator>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container, class OutputIterator>
+    BOOST_MATH_REQUIRES_ITER(BOOST_MATH_OUTPUT_ITER(typename Container::iterator, OutputIterator))
 inline OutputIterator mode(ExecutionPolicy&& exec, Container & v, OutputIterator output)
 {
     return mode(exec, std::begin(v), std::end(v), output);
 }
 
-template<class ForwardIterator, class OutputIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, class OutputIterator>
+    BOOST_MATH_REQUIRES_ITER(BOOST_MATH_OUTPUT_ITER(ForwardIterator, OutputIterator))
 inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output)
 {
     return mode(std::execution::seq, first, last, output);
@@ -628,13 +632,14 @@ inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIt
 
 // Requires enable_if_t to not clash with impl that returns std::list
 // Very ugly. std::is_execution_policy_v returns false for the std::execution objects and decltype of the objects (e.g. std::execution::seq)
-template<class Container, class OutputIterator, std::enable_if_t<!std::is_convertible_v<std::execution::sequenced_policy, Container> &&
+template<BOOST_MATH_CONTAINER Container, class OutputIterator, std::enable_if_t<!std::is_convertible_v<std::execution::sequenced_policy, Container> &&
                                                                  !std::is_convertible_v<std::execution::parallel_unsequenced_policy, Container> &&
                                                                  !std::is_convertible_v<std::execution::parallel_policy, Container>
                                                                  #if __cpp_lib_execution > 201900
                                                                  && !std::is_convertible_v<std::execution::unsequenced_policy, Container>
                                                                  #endif
                                                                  , bool> = true>
+    BOOST_MATH_REQUIRES_ITER(BOOST_MATH_OUTPUT_ITER(typename Container::iterator, OutputIterator))
 inline OutputIterator mode(Container & v, OutputIterator output)
 {
     return mode(std::execution::seq, std::begin(v), std::end(v), output);
@@ -642,27 +647,28 @@ inline OutputIterator mode(Container & v, OutputIterator output)
 
 // std::list is the return type for the proposed STL stats library
 
-template<class ExecutionPolicy, class ForwardIterator, class Real = typename std::iterator_traits<ForwardIterator>::value_type>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_FORWARD_ITER ForwardIterator, 
+         BOOST_MATH_ARBITRARY_NUMERICAL T = typename std::iterator_traits<ForwardIterator>::value_type>
 inline auto mode(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last)
 {
-    std::list<Real> modes;
+    std::list<T> modes;
     mode(exec, first, last, std::inserter(modes, modes.begin()));
     return modes;
 }
 
-template<class ExecutionPolicy, class Container>
+template<BOOST_MATH_EXECUTION_POLICY ExecutionPolicy, BOOST_MATH_CONTAINER Container>
 inline auto mode(ExecutionPolicy&& exec, Container & v)
 {
     return mode(exec, std::begin(v), std::end(v));
 }
 
-template<class ForwardIterator>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator>
 inline auto mode(ForwardIterator first, ForwardIterator last)
 {
     return mode(std::execution::seq, first, last);
 }
 
-template<class Container>
+template<BOOST_MATH_CONTAINER Container>
 inline auto mode(Container & v)
 {
     return mode(std::execution::seq, std::begin(v), std::end(v));
@@ -677,7 +683,7 @@ namespace boost { namespace math { namespace statistics {
 template<bool B, class T = void>
 using enable_if_t = typename std::enable_if<B, T>::type;
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double mean(const ForwardIterator first, const ForwardIterator last)
 {
@@ -685,14 +691,14 @@ inline double mean(const ForwardIterator first, const ForwardIterator last)
     return detail::mean_sequential_impl<double>(first, last);
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double mean(const Container& c)
 {
     return mean(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real mean(const ForwardIterator first, const ForwardIterator last)
 {
@@ -700,28 +706,28 @@ inline Real mean(const ForwardIterator first, const ForwardIterator last)
     return detail::mean_sequential_impl<Real>(first, last);
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real mean(const Container& c)
 {
     return mean(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double variance(const ForwardIterator first, const ForwardIterator last)
 {
     return std::get<2>(detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last));
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double variance(const Container& c)
 {
     return variance(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real variance(const ForwardIterator first, const ForwardIterator last)
 {
@@ -729,14 +735,14 @@ inline Real variance(const ForwardIterator first, const ForwardIterator last)
 
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real variance(const Container& c)
 {
     return variance(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double sample_variance(const ForwardIterator first, const ForwardIterator last)
 {
@@ -745,14 +751,14 @@ inline double sample_variance(const ForwardIterator first, const ForwardIterator
     return n*variance(first, last)/(n-1);
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double sample_variance(const Container& c)
 {
     return sample_variance(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real sample_variance(const ForwardIterator first, const ForwardIterator last)
 {
@@ -761,14 +767,14 @@ inline Real sample_variance(const ForwardIterator first, const ForwardIterator l
     return n*variance(first, last)/(n-1);
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real sample_variance(const Container& c)
 {
     return sample_variance(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline std::pair<double, double> mean_and_sample_variance(const ForwardIterator first, const ForwardIterator last)
 {
@@ -776,14 +782,14 @@ inline std::pair<double, double> mean_and_sample_variance(const ForwardIterator
     return std::make_pair(std::get<0>(results), std::get<3>(results)*std::get<2>(results)/(std::get<3>(results)-1.0));
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline std::pair<double, double> mean_and_sample_variance(const Container& c)
 {
     return mean_and_sample_variance(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline std::pair<Real, Real> mean_and_sample_variance(const ForwardIterator first, const ForwardIterator last)
 {
@@ -791,14 +797,14 @@ inline std::pair<Real, Real> mean_and_sample_variance(const ForwardIterator firs
     return std::make_pair(std::get<0>(results), std::get<3>(results)*std::get<2>(results)/(std::get<3>(results)-Real(1)));
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline std::pair<Real, Real> mean_and_sample_variance(const Container& c)
 {
     return mean_and_sample_variance(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline std::tuple<double, double, double, double> first_four_moments(const ForwardIterator first, const ForwardIterator last)
 {
@@ -807,14 +813,14 @@ inline std::tuple<double, double, double, double> first_four_moments(const Forwa
                            std::get<3>(results) / std::get<4>(results));
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline std::tuple<double, double, double, double> first_four_moments(const Container& c)
 {
     return first_four_moments(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline std::tuple<Real, Real, Real, Real> first_four_moments(const ForwardIterator first, const ForwardIterator last)
 {
@@ -823,42 +829,42 @@ inline std::tuple<Real, Real, Real, Real> first_four_moments(const ForwardIterat
                            std::get<3>(results) / std::get<4>(results));
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline std::tuple<Real, Real, Real, Real> first_four_moments(const Container& c)
 {
     return first_four_moments(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double skewness(const ForwardIterator first, const ForwardIterator last)
 {
     return detail::skewness_sequential_impl<double>(first, last);
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double skewness(const Container& c)
 {
     return skewness(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real skewness(const ForwardIterator first, const ForwardIterator last)
 {
     return detail::skewness_sequential_impl<Real>(first, last);
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real skewness(const Container& c)
 {
     return skewness(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double kurtosis(const ForwardIterator first, const ForwardIterator last)
 {
@@ -874,14 +880,14 @@ inline double kurtosis(const ForwardIterator first, const ForwardIterator last)
     }
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double kurtosis(const Container& c)
 {
     return kurtosis(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real kurtosis(const ForwardIterator first, const ForwardIterator last)
 {
@@ -897,42 +903,42 @@ inline Real kurtosis(const ForwardIterator first, const ForwardIterator last)
     }
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real kurtosis(const Container& c)
 {
     return kurtosis(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double excess_kurtosis(const ForwardIterator first, const ForwardIterator last)
 {
     return kurtosis(first, last) - 3;
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double excess_kurtosis(const Container& c)
 {
     return excess_kurtosis(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real excess_kurtosis(const ForwardIterator first, const ForwardIterator last)
 {
     return kurtosis(first, last) - 3;
 }
 
-template<class Container, typename Real = typename Container::value_type,
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real excess_kurtosis(const Container& c)
 {
     return excess_kurtosis(std::begin(c), std::end(c));
 }
 
-template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type>
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<RandomAccessIterator>::value_type>
 Real median(RandomAccessIterator first, RandomAccessIterator last)
 {
     const auto num_elems = std::distance(first, last);
@@ -952,13 +958,13 @@ Real median(RandomAccessIterator first, RandomAccessIterator last)
     }
 }
 
-template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type>
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename RandomAccessContainer::value_type>
 inline Real median(RandomAccessContainer& c)
 {
     return median(std::begin(c), std::end(c));
 }
 
-template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 {
@@ -970,14 +976,14 @@ inline double gini_coefficient(RandomAccessIterator first, RandomAccessIterator
     return detail::gini_coefficient_sequential_impl<double>(first, last);
 }
 
-template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename RandomAccessContainer::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double gini_coefficient(RandomAccessContainer& c)
 {
     return gini_coefficient(std::begin(c), std::end(c));
 }
 
-template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 {
@@ -989,14 +995,14 @@ inline Real gini_coefficient(RandomAccessIterator first, RandomAccessIterator la
     return detail::gini_coefficient_sequential_impl<Real>(first, last);
 }
 
-template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename RandomAccessContainer::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real gini_coefficient(RandomAccessContainer& c)
 {
     return gini_coefficient(std::begin(c), std::end(c));
 }
 
-template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 {
@@ -1004,14 +1010,14 @@ inline double sample_gini_coefficient(RandomAccessIterator first, RandomAccessIt
     return n*gini_coefficient(first, last)/(n-1);
 }
 
-template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename RandomAccessContainer::value_type,
          enable_if_t<std::is_integral<Real>::value, bool> = true>
 inline double sample_gini_coefficient(RandomAccessContainer& c)
 {
     return sample_gini_coefficient(std::begin(c), std::end(c));
 }
 
-template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<RandomAccessIterator>::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 {
@@ -1019,14 +1025,14 @@ inline Real sample_gini_coefficient(RandomAccessIterator first, RandomAccessIter
     return n*gini_coefficient(first, last)/(n-1);
 }
 
-template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type,
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename RandomAccessContainer::value_type,
          enable_if_t<!std::is_integral<Real>::value, bool> = true>
 inline Real sample_gini_coefficient(RandomAccessContainer& c)
 {
     return sample_gini_coefficient(std::begin(c), std::end(c));
 }
 
-template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type>
+template<BOOST_MATH_RANDOM_ACCESS_ITER RandomAccessIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<RandomAccessIterator>::value_type>
 Real median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last,
     typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())
 {
@@ -1054,14 +1060,14 @@ Real median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator
     }
 }
 
-template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type>
+template<BOOST_MATH_RANDOM_ACCESS_CONTAINER RandomAccessContainer, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename RandomAccessContainer::value_type>
 inline Real median_absolute_deviation(RandomAccessContainer& c,
     typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())
 {
     return median_absolute_deviation(std::begin(c), std::end(c), center);
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type>
 Real interquartile_range(ForwardIterator first, ForwardIterator last)
 {
     static_assert(!std::is_integral<Real>::value, "Integer values have not yet been implemented.");
@@ -1103,13 +1109,13 @@ Real interquartile_range(ForwardIterator first, ForwardIterator last)
     }
 }
 
-template<class Container, typename Real = typename Container::value_type>
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type>
 Real interquartile_range(Container& c)
 {
     return interquartile_range(std::begin(c), std::end(c));
 }
 
-template<class ForwardIterator, class OutputIterator,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, class OutputIterator,
     enable_if_t<std::is_same<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>::value, bool> = true>
 inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output)
 {
@@ -1121,7 +1127,7 @@ inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIt
     return detail::mode_impl(first, last, output);
 }
 
-template<class ForwardIterator, class OutputIterator,
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, class OutputIterator,
     enable_if_t<!std::is_same<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>::value, bool> = true>
 inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output)
 {
@@ -1133,13 +1139,13 @@ inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIt
     return detail::mode_impl(first, last, output);
 }
 
-template<class Container, class OutputIterator>
+template<BOOST_MATH_CONTAINER Container, class OutputIterator>
 inline OutputIterator mode(Container& c, OutputIterator output)
 {
     return mode(std::begin(c), std::end(c), output);
 }
 
-template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type>
+template<BOOST_MATH_FORWARD_ITER ForwardIterator, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename std::iterator_traits<ForwardIterator>::value_type>
 inline std::list<Real> mode(ForwardIterator first, ForwardIterator last)
 {
     std::list<Real> modes;
@@ -1147,7 +1153,7 @@ inline std::list<Real> mode(ForwardIterator first, ForwardIterator last)
     return modes;
 }
 
-template<class Container, typename Real = typename Container::value_type>
+template<BOOST_MATH_CONTAINER Container, BOOST_MATH_ARBITRARY_NUMERICAL Real = typename Container::value_type>
 inline std::list<Real> mode(Container& c)
 {
     return mode(std::begin(c), std::end(c));
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 9ab16b3da3..68d631c9fb 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1647,6 +1647,46 @@ alias no_eh_tests :
 
 explicit no_eh_tests ;
 
+# alias for the ccmath tests
+
+alias ccmath_tests :
+  ccmath_abs_test
+  ccmath_ceil_test
+  ccmath_copysign_test
+  ccmath_div_test
+  ccmath_fdim_test
+  ccmath_floor_test
+  ccmath_fma_test
+  ccmath_fmax_test
+  ccmath_fmin_test
+  ccmath_fmod_test
+  ccmath_fpclassify_test
+  ccmath_frexp_test
+  ccmath_hypot_test
+  ccmath_ilogb_test
+  ccmath_isfinite_test
+  ccmath_isgreater_test
+  ccmath_isgreaterequal_test
+  ccmath_isinf_test
+  ccmath_isless_test
+  ccmath_islessequal_test
+  ccmath_isnan_test
+  ccmath_isnormal_test
+  ccmath_isunordered_test
+  ccmath_ldexp_test
+  ccmath_logb_test
+  ccmath_modf_test
+  ccmath_next_test
+  ccmath_remainder_test
+  ccmath_round_test
+  ccmath_scalbln_test
+  ccmath_scalbn_test
+  ccmath_signbit_test
+  ccmath_sqrt_test
+;
+
+explicit ccmath_tests ;
+
 # Some aliases which group blocks of tests for CI testing:
 
 alias github_ci_block_1 : special_fun float128_tests distribution_tests mp misc ;