diff --git a/include/boost/math/distributions/binomial.hpp b/include/boost/math/distributions/binomial.hpp
index 48f5f016f9..92e24b96bf 100644
--- a/include/boost/math/distributions/binomial.hpp
+++ b/include/boost/math/distributions/binomial.hpp
@@ -822,8 +822,8 @@ BOOST_MATH_NAMESPACE_BEGIN
       // denominator is (n + 1) p - (k + 1):
       RealType k_right = mode;
       RealType k_left = mode;
-      detail::entropy_compensated_sum<RealType> s_right;
-      detail::entropy_compensated_sum<RealType> s_left;
+      detail::compensated_sum<RealType> s_right;
+      detail::compensated_sum<RealType> s_left;
       auto log_ratio = [&](RealType k, RealType& ratio, bool left)
       {
         return detail::entropy_log_ratio(RealType((n - k) * p), RealType((k + 1) * q), RealType(fma(n + 1, p, -(k + 1))), ratio, left, Policy());
diff --git a/include/boost/math/distributions/detail/entropy.hpp b/include/boost/math/distributions/detail/entropy.hpp
index 77d019eb17..660d375108 100644
--- a/include/boost/math/distributions/detail/entropy.hpp
+++ b/include/boost/math/distributions/detail/entropy.hpp
@@ -16,7 +16,7 @@
 #include <boost/math/special_functions/digamma.hpp>
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/special_functions/log1p.hpp>
-#include <boost/math/special_functions/detail/two_product.hpp>
+#include <boost/math/tools/compensated_arithmetic.hpp>
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/policies/error_handling.hpp>
 #ifndef BOOST_MATH_BUILD_MODULE
@@ -126,34 +126,6 @@ inline RealType inverse_gamma_entropy_unit_scale(RealType a, const Policy& pol)
    return constants::log_root_two_pi<RealType>() + constants::half<RealType>() - 3 * log(a) / 2 + 1 / (2 * a) + r + 2 * nu;
 }
 
-// a * b - c * d for a * b close to c * d: with p1 + e1 = a * b and p2 + e2 = c * d exactly (two_product,
-// which also works for types whose fma isn't fused), p1 - p2 is exact by Sterbenz's lemma.
-template <class RealType>
-inline RealType entropy_difference_of_products(RealType a, RealType b, RealType c, RealType d)
-{
-   two_product<RealType> product;
-   RealType p1, e1, p2, e2;
-   product(a, b, p1, e1);
-   product(c, d, p2, e2);
-   return (p1 - p2) + (e1 - e2);
-}
-
-// Neumaier's compensated sum: the direct sums can run to ~1e5 terms (negative binomials with small p).
-template <class RealType>
-struct entropy_compensated_sum
-{
-   RealType sum = 0;
-   RealType c = 0;
-   void add(RealType x)
-   {
-      BOOST_MATH_STD_USING
-      RealType t = sum + x;
-      c += (fabs(sum) >= fabs(x)) ? RealType((sum - t) + x) : RealType((x - t) + sum);
-      sum = t;
-   }
-   RealType value() const { return sum + c; }
-};
-
 // Shannon entropy -sum p log(p) of a discrete distribution, from the log ratios S = log(p / p_0) of its
 // probabilities at points spaced h apart, relative to a point near the mode (S = 0 there):
 //    H = log(h Z) - sum(w S) / Z,  where w = exp(S) and Z = sum(w).
@@ -175,8 +147,9 @@ RealType discrete_entropy_sum(Right& right, Left& left, RealType h, RealType rho
    BOOST_MATH_STD_USING
    const RealType eps = policies::get_epsilon<RealType, Policy>();
    const std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
-   entropy_compensated_sum<RealType> z;  // sum of w, without the w = 1 at S = 0.
-   entropy_compensated_sum<RealType> zs; // sum of w S.
+   // Neumaier's compensated sums: the direct sums can run to ~1e5 terms (negative binomials with small p).
+   compensated_sum<RealType> z;  // sum of w, without the w = 1 at S = 0.
+   compensated_sum<RealType> zs; // sum of w S.
    std::uintmax_t iter = 0;
    for (int side = 0; side < 2; ++side)
    {
diff --git a/include/boost/math/distributions/hypergeometric.hpp b/include/boost/math/distributions/hypergeometric.hpp
index 348856d521..6d51d1aa4b 100644
--- a/include/boost/math/distributions/hypergeometric.hpp
+++ b/include/boost/math/distributions/hypergeometric.hpp
@@ -412,7 +412,7 @@ BOOST_MATH_NAMESPACE_BEGIN
       if ((sigma >= 2 * sigmas) && ((std::min)((std::min)(z1, z2), (std::min)(z3, z4)) >= margin))
       {
          // Trapezoidal rule on the continuation, with linear coefficient alpha = log(z2 z3 / (z1 z4)):
-         RealType alpha = BOOST_MATH_NAMESPACE::log1p(RealType(detail::entropy_difference_of_products(z2, z3, z1, z4) / (z1 * z4)), Policy());
+         RealType alpha = BOOST_MATH_NAMESPACE::log1p(RealType(detail::difference_of_products(z2, z3, z1, z4) / (z1 * z4)), Policy());
          RealType h = sigma / detail::discrete_entropy_step_ratio<RealType, Policy>();
          auto s_of = [&](RealType t)
          {
@@ -443,14 +443,14 @@ BOOST_MATH_NAMESPACE_BEGIN
       // whose numerator minus denominator is (r + 1)(n + 1) - (k + 1)(N + 2):
       std::uint64_t k_right = mode;
       std::uint64_t k_left = mode;
-      detail::entropy_compensated_sum<RealType> s_right;
-      detail::entropy_compensated_sum<RealType> s_left;
+      detail::compensated_sum<RealType> s_right;
+      detail::compensated_sum<RealType> s_left;
       auto log_ratio = [&](std::uint64_t k, RealType& ratio, bool left)
       {
          RealType kk = static_cast<RealType>(k);
          RealType num = (rr - kk) * (nn - kk);
          RealType den = (kk + 1) * static_cast<RealType>(N + k + 1 - r - n);
-         RealType d = detail::entropy_difference_of_products(RealType(rr + 1), RealType(nn + 1), RealType(kk + 1), RealType(NN + 2));
+         RealType d = detail::difference_of_products(RealType(rr + 1), RealType(nn + 1), RealType(kk + 1), RealType(NN + 2));
          return detail::entropy_log_ratio(num, den, d, ratio, left, Policy());
       };
       auto right = [&](RealType& s, RealType& ratio)
diff --git a/include/boost/math/distributions/kolmogorov_smirnov.hpp b/include/boost/math/distributions/kolmogorov_smirnov.hpp
index 891f0cafbf..c4ae33ff8d 100644
--- a/include/boost/math/distributions/kolmogorov_smirnov.hpp
+++ b/include/boost/math/distributions/kolmogorov_smirnov.hpp
@@ -123,6 +123,7 @@
 #include <boost/math/distributions/detail/common_error_handling.hpp>
 #include <boost/math/special_functions/jacobi_theta.hpp>
 #include <boost/math/special_functions/log1p.hpp>
+#include <boost/math/tools/compensated_arithmetic.hpp>
 #include <boost/math/tools/big_constant.hpp>
 #include <boost/math/tools/tuple.hpp>
 #include <boost/math/tools/roots.hpp> // Newton-Raphson
@@ -160,42 +161,16 @@ inline RealType kolmogorov_smirnov_kurtosis_excess_constant() {
     return BOOST_MATH_BIG_CONSTANT(RealType, 1000, 0.8816189679105236704015538304328295408355055629289034536289160037385510119639244078254098968815665536);
 }
 
-// Splits t = hi + lo so that hi carries only the leading half of the
-// significand: products of two such hi parts are then exact (Veltkamp's
-// splitting, with splitter = 2^ceil(digits/2) + 1).
-template <class RealType>
-inline RealType kolmogorov_smirnov_splitter() {
-    BOOST_MATH_STD_USING
-    return ldexp(RealType(1), (tools::digits<RealType>() + 1) / 2) + 1;
-}
-
-template <class RealType>
-inline void kolmogorov_smirnov_split(RealType t, RealType splitter, RealType& hi, RealType& lo) {
-    RealType g = splitter * t;
-    hi = g - (g - t);
-    lo = t - hi;
-}
-
-// Returns the rounding error of the product ab = a*b, i.e. the exact product
-// is ab + result (Dekker's algorithm, which works for any floating-point
-// type, unlike fma).
-template <class RealType>
-inline RealType kolmogorov_smirnov_product_error(RealType a, RealType b, RealType ab, RealType splitter) {
-    RealType ah, al, bh, bl;
-    kolmogorov_smirnov_split(a, splitter, ah, al);
-    kolmogorov_smirnov_split(b, splitter, bh, bl);
-    return (((ah * bh - ab) + ah * bl) + al * bh) + al * bl;
-}
-
 // Given u = x*x*n as computed in floating point, returns the rounding error
 // in u, so that exp(-c*u) can be evaluated as exp(-c*u) * (1 - c*err)
 // without the rounding error of u being amplified by the size of the
 // exponent c*u.
 template <class RealType>
-inline RealType kolmogorov_smirnov_x2n_error(RealType x, RealType n, RealType u, RealType splitter) {
+inline RealType kolmogorov_smirnov_x2n_error(RealType x, RealType n, RealType u) {
+    two_product<RealType> product;
     RealType x2 = x * x;
-    RealType e1 = kolmogorov_smirnov_product_error(x, x, x2, splitter);
-    RealType e2 = kolmogorov_smirnov_product_error(x2, n, u, splitter);
+    RealType e1 = product.error(x, x, x2);
+    RealType e2 = product.error(x2, n, u);
     return e2 + e1 * n; // x*x*n == u + err
 }
 
@@ -211,7 +186,7 @@ inline RealType kolmogorov_smirnov_nome(RealType x, RealType n) {
     RealType q = exp(-2 * u);
     if (q == 0)
         return q;
-    RealType err = kolmogorov_smirnov_x2n_error(x, n, u, kolmogorov_smirnov_splitter<RealType>());
+    RealType err = kolmogorov_smirnov_x2n_error(x, n, u);
     return q * (1 - 2 * err); // exp(-2 (u + err))
 }
 
@@ -262,10 +237,10 @@ RealType kolmogorov_smirnov_pdf_small_x(RealType x, RealType n, const Policy&) {
     // The rounding errors in the constant, in u and in the division, to
     // second order, so that exp(-W) can be corrected by a factor (1 - delta):
     //   W_exact = W + delta,  delta = ((c_hi + c_lo) - W*u - W*err) / u,  u_exact = u + err
-    RealType splitter = kolmogorov_smirnov_splitter<RealType>();
-    RealType err = kolmogorov_smirnov_x2n_error(x, n, u, splitter);
-    RealType Wu = W * u;
-    RealType delta = ((((c_hi - Wu) - kolmogorov_smirnov_product_error(W, u, Wu, splitter)) + c_lo) - W * err) / u;
+    RealType err = kolmogorov_smirnov_x2n_error(x, n, u);
+    RealType Wu, eWu;
+    two_product<RealType>()(W, u, Wu, eWu);
+    RealType delta = ((((c_hi - Wu) - eWu) + c_lo) - W * err) / u;
     if (half != 0) {
         RealType leading = half * (1 - delta / 2) * sqrt(constants::root_two_pi<RealType>() * sqrt(n) * (pi2 / 4 - u)) / u;
         return leading * leading;
@@ -305,11 +280,11 @@ inline RealType kolmogorov_smirnov_pdf_large_x(RealType x, RealType n, const Pol
         RealType half = exp(-u);
         if (half == 0)
             return static_cast<RealType>(0);
-        RealType err = kolmogorov_smirnov_x2n_error(x, n, u, kolmogorov_smirnov_splitter<RealType>());
+        RealType err = kolmogorov_smirnov_x2n_error(x, n, u);
         RealType leading = half * (1 - err) * sqrt(8 * x * n);
         return leading * leading;
     }
-    RealType err = kolmogorov_smirnov_x2n_error(x, n, u, kolmogorov_smirnov_splitter<RealType>());
+    RealType err = kolmogorov_smirnov_x2n_error(x, n, u);
     r *= (1 - 2 * err); // exp(-2 (u + err))
     RealType r2 = r * r;
 
diff --git a/include/boost/math/distributions/negative_binomial.hpp b/include/boost/math/distributions/negative_binomial.hpp
index a700629800..f12d966a9e 100644
--- a/include/boost/math/distributions/negative_binomial.hpp
+++ b/include/boost/math/distributions/negative_binomial.hpp
@@ -654,8 +654,7 @@ BOOST_MATH_NAMESPACE_BEGIN
         // z1 q - z2 = (r - 1) - z1 p. z1 = z1_hi + z1_lo exactly; z1_lo enters to first order.
         RealType d = r - 1;
         RealType z1_hi = z2 + d;
-        RealType bb = z1_hi - z2;
-        RealType z1_lo = (z2 - (z1_hi - bb)) + (d - bb);
+        RealType z1_lo = detail::two_sum_error(z2, d, z1_hi);
         RealType alpha = BOOST_MATH_NAMESPACE::log1p(RealType((fma(-z1_hi, p, d) - z1_lo * p) / z2), Policy());
         RealType h = sigma / detail::discrete_entropy_step_ratio<RealType, Policy>();
         auto s_of = [&](RealType t)
@@ -680,8 +679,8 @@ BOOST_MATH_NAMESPACE_BEGIN
       // denominator is r q - 1 - k p. For r < 1 the ratio rises towards q in the tail:
       RealType k_right = mode;
       RealType k_left = mode;
-      detail::entropy_compensated_sum<RealType> s_right;
-      detail::entropy_compensated_sum<RealType> s_left;
+      detail::compensated_sum<RealType> s_right;
+      detail::compensated_sum<RealType> s_left;
       auto log_ratio = [&](RealType k, RealType& ratio, bool left)
       {
         return detail::entropy_log_ratio(RealType((k + r) * q), RealType(k + 1), RealType(fma(r, q, -fma(k, p, RealType(1)))), ratio, left, Policy());
diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp
index c90cb51824..b20d48ccea 100644
--- a/include/boost/math/distributions/poisson.hpp
+++ b/include/boost/math/distributions/poisson.hpp
@@ -650,8 +650,8 @@ BOOST_MATH_NAMESPACE_BEGIN
       RealType mode = floor(lambda);
       RealType k_right = mode;
       RealType k_left = mode;
-      detail::entropy_compensated_sum<RealType> s_right;
-      detail::entropy_compensated_sum<RealType> s_left;
+      detail::compensated_sum<RealType> s_right;
+      detail::compensated_sum<RealType> s_left;
       auto right = [&](RealType& s, RealType& ratio)
       {
         s_right.add(detail::entropy_log_ratio(lambda, RealType(k_right + 1), RealType(lambda - (k_right + 1)), ratio, false, Policy()));
diff --git a/include/boost/math/quadrature/gauss_hermite.hpp b/include/boost/math/quadrature/gauss_hermite.hpp
index 7b593ed91e..ec2d7ef2f4 100644
--- a/include/boost/math/quadrature/gauss_hermite.hpp
+++ b/include/boost/math/quadrature/gauss_hermite.hpp
@@ -30,6 +30,7 @@
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/policies/policy.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/tools/compensated_arithmetic.hpp>
 #include <boost/math/tools/big_constant.hpp>
 #include <boost/math/tools/precision.hpp>
 #include <boost/math/quadrature/detail/quadrature_constant.hpp>
@@ -189,16 +190,12 @@ class gauss_hermite_detail
       // their weighted mean, plus two roundings. That removes the shared part, which is what integrals of
       // smooth functions see, and moves no weight by more than |e| <= max |e_i|; measured, both the largest and
       // the mean |e_i| fall for N = 200, 500 and 706. Sum from the smallest weight, with Neumaier's compensation:
-      Real sum = 0;
-      Real compensation = 0;
+      BOOST_MATH_NAMESPACE::detail::compensated_sum<Real> sum;
       for (std::size_t i = w.size(); i-- > 0;)
       {
-         Real term = (offset && (i == 0)) ? w[i] : Real(2 * w[i]);
-         Real t = sum + term;
-         compensation += abs(sum) >= abs(term) ? Real((sum - t) + term) : Real((term - t) + sum);
-         sum = t;
+         sum.add((offset && (i == 0)) ? w[i] : Real(2 * w[i]));
       }
-      const Real scale = BOOST_MATH_NAMESPACE::constants::root_pi<Real>() / (sum + compensation);
+      const Real scale = BOOST_MATH_NAMESPACE::constants::root_pi<Real>() / sum.value();
       for (Real& weight_i : w)
       {
          weight_i *= scale;
diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index dfe19ee4b4..998b0a6495 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/include/boost/math/special_functions/beta.hpp
@@ -31,9 +31,7 @@
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/binomial.hpp>
 #include <boost/math/special_functions/factorials.hpp>
-#ifndef BOOST_MATH_HAS_GPU_SUPPORT
-#include <boost/math/special_functions/detail/two_product.hpp>  // only the host-only no-Lanczos power terms use it
-#endif
+#include <boost/math/tools/compensated_arithmetic.hpp>
 #include <boost/math/tools/roots.hpp>
 
 BOOST_MATH_NAMESPACE_BEGIN
@@ -369,16 +367,6 @@ BOOST_MATH_GPU_ENABLED T lbeta_imp(T a, T b, const lanczos::undefined_lanczos& l
 } // template <class T>T lbeta_imp(T a, T b, const lanczos::undefined_lanczos& l)
 #endif
 
-//
-// The exact sum u + v is s + (return value), where s = fl(u + v):
-//
-template <class T>
-BOOST_MATH_GPU_ENABLED inline T ibeta_sum_error(T u, T v, T s)
-{
-   T vv = s - u;
-   return (u - (s - vv)) + (v - vv);
-}
-
 //
 // The bases of the power terms minus one, l1 = x cgh / agh - 1 and l2 = y cgh / bgh - 1, where
 // agh = a + gh, bgh = b + gh and cgh = a + b + gh (gh = 0 for no shift).
@@ -395,9 +383,9 @@ BOOST_MATH_GPU_ENABLED inline void ibeta_power_bases(T a, T b, T x, T y, T gh, T
    using std::fma;
    #endif
    T c = a + b;
-   T agh_err = ibeta_sum_error(a, gh, agh);
-   T bgh_err = ibeta_sum_error(b, gh, bgh);
-   T cgh_err = ibeta_sum_error(c, gh, cgh) + ibeta_sum_error(a, b, c);
+   T agh_err = two_sum_error(a, gh, agh);
+   T bgh_err = two_sum_error(b, gh, bgh);
+   T cgh_err = two_sum_error(c, gh, cgh) + two_sum_error(a, b, c);
    if(x <= y)
    {
       l1 = (fma(x, cgh, -agh) + (x * cgh_err - agh_err)) / agh;
@@ -422,7 +410,7 @@ BOOST_MATH_GPU_ENABLED inline T ibeta_power_log_near_mean(T a, T b, T l1, T l2,
    #ifndef BOOST_MATH_HAS_GPU_SUPPORT
    using std::fma;
    #endif
-   T bgh_err = ibeta_sum_error(b, gh, bgh);
+   T bgh_err = two_sum_error(b, gh, bgh);
    //
    // Both bases are near 1, and when a and b are large the two terms of the exponent
    // a log1p(l1) + b log1p(l2) nearly cancel. Since agh (1 + l1) + bgh (1 + l2) = cgh,
@@ -435,11 +423,11 @@ BOOST_MATH_GPU_ENABLED inline T ibeta_power_log_near_mean(T a, T b, T l1, T l2,
    // rounding error found with fma, is lh + ll, and the rounding of lh + the rest is added back below.
    //
    T ab = a - b;
-   T ab_err = ibeta_sum_error(a, T(-b), ab);
+   T ab_err = two_sum_error(a, T(-b), ab);
    T ph = fma(l1, ab, T(0));  // fma(x, y, 0) so that the product is not contracted into a later sum
    T pl = fma(l1, ab, -ph) + l1 * ab_err;
    T nh = ph - b;
-   T nl = pl + ibeta_sum_error(ph, T(-b), nh);
+   T nl = pl + two_sum_error(ph, T(-b), nh);
    T wh = nh / bgh;
    T wl = (fma(-wh, bgh, nh) + (nl - wh * bgh_err)) / bgh;
    T lh = fma(gh, wh, T(0));
@@ -448,7 +436,7 @@ BOOST_MATH_GPU_ENABLED inline T ibeta_power_log_near_mean(T a, T b, T l1, T l2,
       + a * BOOST_MATH_NAMESPACE::log1pmx(l1, pol)
       + b * BOOST_MATH_NAMESPACE::log1pmx(l2, pol);
    T l = lh + rest;
-   l_err = ibeta_sum_error(lh, rest, l);
+   l_err = two_sum_error(lh, rest, l);
    return l;
 }
 
@@ -731,7 +719,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a,
       // from the smaller of x and y and the rounding error of c, the bases are 1 + d / a and 1 - d / b,
       // and the d terms of the exponent cancel exactly, leaving a log1pmx(d / a) + b log1pmx(-d / b):
       //
-      T c_err = ibeta_sum_error(a, b, c);
+      T c_err = two_sum_error(a, b, c);
       T ph, pl;
       two_product<T> product;
       product(x <= y ? x : y, c, ph, pl);  // ph is fma(x, c, 0), which contraction can't fuse into ph - a
@@ -879,22 +867,16 @@ BOOST_MATH_GPU_ENABLED T ibeta_series_sum(ibeta_series_t<T>& s, T s0, BOOST_MATH
    BOOST_MATH_STD_USING
    const T factor = BOOST_MATH_NAMESPACE::policies::get_epsilon<T, Policy>();
    BOOST_MATH_NAMESPACE::uintmax_t counter = max_iter;
-   T sum = s0;
-   T comp = 0;
+   compensated_sum<T> sum(s0);
    T term;
    do
    {
       term = s();
-      T next = sum + term;
-      if (abs(sum) >= abs(term))
-         comp += (sum - next) + term;
-      else
-         comp += (term - next) + sum;
-      sum = next;
+      sum.add(term);
    }
-   while ((abs(factor * sum) < abs(term)) && --counter);
+   while ((abs(factor * sum.sum) < abs(term)) && --counter);
    max_iter = max_iter - counter;
-   return sum + comp;
+   return sum.value();
 }
 
 template <class T, class Lanczos, class Policy>
@@ -1297,9 +1279,8 @@ BOOST_MATH_GPU_ENABLED T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T
    //
    // The sum is often started from s0 = -1, when the first term is of order one and each later term
    // is rounded onto a coarse grid, so keep their rounding errors in comp (Neumaier summation):
-   T pj = prefix * j;
-   T sum = s0 + pj;  // Value at N = 0
-   T comp = (fabs(s0) >= fabs(pj)) ? T((s0 - sum) + pj) : T((pj - sum) + s0);
+   compensated_sum<T> sum(s0);
+   sum.add(prefix * j);  // Value at N = 0
    // some variables we'll need:
    unsigned tnp1 = 1; // 2*N+1
    T lx2 = lx / 2;
@@ -1346,15 +1327,13 @@ BOOST_MATH_GPU_ENABLED T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T
       // pull it together with Eq 9:
       //
       T r = prefix * p[n] * j;
-      T next = sum + r;
-      comp += (fabs(sum) >= fabs(r)) ? T((sum - next) + r) : T((r - next) + sum);
-      sum = next;
+      sum.add(r);
       // r is always small:
       BOOST_MATH_ASSERT(tools::max_value<T>() * tools::epsilon<T>() > fabs(r));
-      if(fabs(r / tools::epsilon<T>()) < fabs(sum))
+      if(fabs(r / tools::epsilon<T>()) < fabs(sum.sum))
          break;
    }
-   return sum + comp;
+   return sum.value();
 } // template <class T, class Lanczos>T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Lanczos& l, bool normalised)
 
 //
diff --git a/include/boost/math/special_functions/detail/two_product.hpp b/include/boost/math/special_functions/detail/two_product.hpp
deleted file mode 100644
index ae6752919f..0000000000
--- a/include/boost/math/special_functions/detail/two_product.hpp
+++ /dev/null
@@ -1,99 +0,0 @@
-//  (C) Copyright Nick Thompson 2026.
-//  Use, modification and distribution are subject to the
-//  Boost Software License, Version 1.0. (See accompanying file
-//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_TWO_PRODUCT_HPP
-#define BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_TWO_PRODUCT_HPP
-
-#ifdef _MSC_VER
-#pragma once
-#endif
-
-#include <boost/math/tools/config.hpp>
-#include <boost/math/tools/precision.hpp>
-
-#ifndef BOOST_MATH_BUILD_MODULE
-#include <cmath>
-#endif
-
-BOOST_MATH_NAMESPACE_BEGIN namespace detail {
-
-namespace two_product_fma {
-using std::fma;
-
-// result = fma(a, b, c) where the type has an fma, found by argument-dependent lookup for the
-// non-built-in types.  Returns false, leaving result alone, where there is none.
-template <class T>
-inline auto evaluate(const T& a, const T& b, const T& c, T& result, int) -> decltype(static_cast<void>(fma(a, b, c)), bool())
-{
-   result = fma(a, b, c);
-   return true;
-}
-template <class T>
-inline bool evaluate(const T&, const T&, const T&, T&, long)
-{
-   return false;
-}
-} // namespace two_product_fma
-
-//
-// Splits a product into its rounded value and its rounding error: p + e = a b exactly.  A fused fma gives e
-// directly, but having an fma is no guarantee of one: the Boost.Multiprecision backends without a native one
-// (cpp_bin_float among them) evaluate it as a * b + c with two roundings, which would make every e zero.  So
-// probe once per type: (1 + eps)^2 = 1 + 2 eps + eps^2, and only a fused fma(x, x, -(1 + 2 eps)) returns the
-// last term.  Otherwise split the factors in half (Veltkamp) so that the partial products are exact (Dekker).
-// The rounded product comes from fma(a, b, 0), which the compiler can't contract into a later addition.
-//
-template <class T>
-class two_product
-{
-public:
-   two_product() : fused(is_fused()), splitter(0)
-   {
-      BOOST_MATH_STD_USING
-      if (!fused)
-         splitter = ldexp(T(1), (tools::digits<T>() + 1) / 2) + 1;
-   }
-   void operator()(const T& a, const T& b, T& p, T& e) const
-   {
-      if (fused)
-      {
-         two_product_fma::evaluate(a, b, T(0), p, 0);
-         two_product_fma::evaluate(a, b, T(-p), e, 0);
-         return;
-      }
-      p = a * b;
-      T ah, al, bh, bl;
-      split(a, ah, al);
-      split(b, bh, bl);
-      e = (((ah * bh - p) + ah * bl) + al * bh) + al * bl;
-   }
-private:
-   void split(const T& t, T& hi, T& lo) const
-   {
-      T g = splitter * t;
-      hi = g - (g - t);
-      lo = t - hi;
-   }
-   static bool is_fused()
-   {
-      static const bool result = []
-      {
-         // From a volatile, so that the compiler can't fold the probe at compile time:
-         volatile int volatile_one = 1;
-         T one = static_cast<int>(volatile_one);
-         T eps = tools::epsilon<T>();
-         T x = one + eps;
-         T eps_squared = 0;
-         return two_product_fma::evaluate(x, x, T(-(one + 2 * eps)), eps_squared, 0) && (eps_squared != 0);
-      }();
-      return result;
-   }
-   bool fused;
-   T splitter;
-};
-
-} BOOST_MATH_NAMESPACE_END // namespace detail
-
-#endif // BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_TWO_PRODUCT_HPP
diff --git a/include/boost/math/special_functions/hypergeometric_2F0.hpp b/include/boost/math/special_functions/hypergeometric_2F0.hpp
index e95350a0ba..5ed64e55b2 100644
--- a/include/boost/math/special_functions/hypergeometric_2F0.hpp
+++ b/include/boost/math/special_functions/hypergeometric_2F0.hpp
@@ -16,7 +16,7 @@
 #include <boost/math/special_functions/laguerre.hpp>
 #include <boost/math/special_functions/hermite.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
-#include <boost/math/special_functions/detail/two_product.hpp>
+#include <boost/math/tools/compensated_arithmetic.hpp>
 #include <boost/math/tools/fraction.hpp>
 
 BOOST_MATH_NAMESPACE_BEGIN namespace detail {
@@ -70,27 +70,23 @@ BOOST_MATH_NAMESPACE_BEGIN namespace detail {
    {
       BOOST_MATH_STD_USING
       two_product<T> product;
-      T t;
       // g1 + e1 = 1 - b z exactly:
       T bz, rbz;
       product(b, z, bz, rbz);
       T g1 = 1 - bz;
-      t = g1 - 1;
-      T e1 = ((1 - (g1 - t)) + (-bz - t)) - rbz;
+      T e1 = two_sum_error(T(1), T(-bz), g1) - rbz;
       T g0 = 1;
       T e0 = 0;
       for (unsigned k = 1; k < n; ++k)
       {
          // a + ea = 1 - (k + b) z and c + ec = k z exactly:
          T u = k + b;
-         t = u - k;
-         T eu = (k - (u - t)) + (b - t);
+         T eu = two_sum_error(T(k), b, u);
          T v, rv;
          product(u, z, v, rv);
          rv += eu * z;
          T a = 1 - v;
-         t = a - 1;
-         T ea = ((1 - (a - t)) + (-v - t)) - rv;
+         T ea = two_sum_error(T(1), T(-v), a) - rv;
          T c, ec;
          product(T(k), z, c, ec);
          // s1 + r1 = a g1, s2 + r2 = c g0 and d + ed = s1 + s2 exactly:
@@ -98,8 +94,7 @@ BOOST_MATH_NAMESPACE_BEGIN namespace detail {
          product(a, g1, s1, r1);
          product(c, g0, s2, r2);
          T d = s1 + s2;
-         t = d - s1;
-         T ed = (s1 - (d - t)) + (s2 - t);
+         T ed = two_sum_error(s1, s2, d);
          T e2 = a * e1 + c * e0 + ed + r1 + r2 + ea * g1 + ec * g0;
          g0 = g1;
          g1 = d;
diff --git a/include/boost/math/special_functions/jacobi_theta.hpp b/include/boost/math/special_functions/jacobi_theta.hpp
index d6cf9c6263..09abfb827f 100644
--- a/include/boost/math/special_functions/jacobi_theta.hpp
+++ b/include/boost/math/special_functions/jacobi_theta.hpp
@@ -109,6 +109,7 @@
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/special_functions/expm1.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/tools/compensated_arithmetic.hpp>
 
 BOOST_MATH_NAMESPACE_BEGIN
 
@@ -237,86 +238,9 @@ struct jacobi_theta_q_power {
 // tau < 1, the square of z plus a multiple of pi/2). Rounding E to working
 // precision costs |E| ulps in the exponential, so instead the rounding
 // errors of pi, of the products and of the division are tracked exactly with
-// Dekker's product and Knuth's sum (which work for any binary floating-point
+// two_product and two_sum_error (which work for any binary floating-point
 // type) and applied as a correction factor exp(-dE) = 1 - dE. What remains
 // is the rounding of tau itself, which is the caller's.
-// (Cached per type: for a type whose precision can change at run time the
-// cached value may split at the wrong place, in which case the corrections
-// below are merely inexact, i.e. no worse than not applying them.)
-template <class RealType>
-inline RealType jacobi_theta_splitter() {
-    BOOST_MATH_STD_USING
-    static const RealType splitter = ldexp(RealType(1), (tools::digits<RealType>() + 1) / 2) + 1;
-    return splitter;
-}
-
-// Veltkamp's splitting: t = hi + lo where hi holds only the leading half of
-// the significand, so that products of hi parts are exact.
-template <class RealType>
-inline void jacobi_theta_split(RealType t, RealType splitter, RealType& hi, RealType& lo) {
-    RealType g = splitter * t;
-    hi = g - (g - t);
-    lo = t - hi;
-}
-
-// result = fma(a, b, c) if the type has an fma, which is std::fma for the
-// built-in types and is found by argument-dependent lookup for the others
-// (the Boost.Multiprecision types all have one). Returns false, leaving
-// result alone, if there is none.
-namespace jacobi_theta_fma {
-using std::fma;
-
-template <class RealType>
-inline auto evaluate(const RealType& a, const RealType& b, const RealType& c, RealType& result, int) -> decltype(static_cast<void>(fma(a, b, c)), bool()) {
-    result = fma(a, b, c);
-    return true;
-}
-template <class RealType>
-inline bool evaluate(const RealType&, const RealType&, const RealType&, RealType&, long) {
-    return false;
-}
-} // namespace jacobi_theta_fma
-
-// Whether the type's fma is actually fused. Having one is no guarantee: some
-// C libraries (newlib, and hence Cygwin) implement fma(x, y, z) as x * y + z
-// with two roundings, as do the Boost.Multiprecision backends without a
-// native one (cpp_bin_float, but not mpfr_float or float128), and that would
-// silently zero every correction below. Probed once per type: (1 + eps)^2 is
-// 1 + 2 eps + eps^2, and the last term is exactly what a fused
-// fma(x, x, -(1 + 2 eps)) returns, while x * x + z rounds it away. The inputs
-// are built from a volatile so that the compiler cannot fold the call at
-// compile time.
-template <class RealType>
-inline bool jacobi_theta_fma_is_fused() {
-    static const bool fused = [] {
-        volatile int volatile_one = 1;
-        RealType one = static_cast<int>(volatile_one), eps = tools::epsilon<RealType>();
-        RealType x = one + eps, z = -(one + 2 * eps), eps_squared = 0;
-        return jacobi_theta_fma::evaluate(x, x, z, eps_squared, 0) && eps_squared != 0;
-    }();
-    return fused;
-}
-
-// The exact product a*b is ab + (return value), where ab = fl(a*b). A fused
-// fma(a, b, -ab) is exactly this quantity, so it is used where the type has
-// one; the splitting is for the types and libraries that do not.
-template <class RealType>
-inline RealType jacobi_theta_product_error(RealType a, RealType b, RealType ab, RealType splitter) {
-    RealType error;
-    if (jacobi_theta_fma_is_fused<RealType>() && jacobi_theta_fma::evaluate(a, b, RealType(-ab), error, 0))
-        return error;
-    RealType ah, al, bh, bl;
-    jacobi_theta_split(a, splitter, ah, al);
-    jacobi_theta_split(b, splitter, bh, bl);
-    return (((ah * bh - ab) + ah * bl) + al * bh) + al * bl;
-}
-
-// The exact sum a+b is s + (return value), where s = fl(a+b).
-template <class RealType>
-inline RealType jacobi_theta_sum_error(RealType a, RealType b, RealType s) {
-    RealType bb = s - a;
-    return (a - (s - bb)) + (b - bb);
-}
 
 // Pi as hi + lo + lo2: hi holds the leading 24 bits exactly, so that hi * e
 // is exact for small integers e, and lo the following ones to the type's own
@@ -358,7 +282,7 @@ inline RealType jacobi_theta_pi(RealType& lo, RealType& lo2) {
 template <class RealType>
 struct jacobi_theta_exponents {
     RealType tau;           // only used for the scale factor 1/sqrt(tau)
-    RealType splitter;
+    two_product<RealType> product;
     RealType pi_hi, pi_lo, pi_lo2;  // pi = pi_hi + pi_lo + pi_lo2
     RealType a_hi, a_lo;    // pi * tau = -ln(q) = a_hi + a_lo
     RealType inv_a_hi;      // 1 / a_hi, only needed when tau < 1
@@ -383,21 +307,23 @@ struct jacobi_theta_exponents {
 
     // pi * e as hi + lo for a small exact e
     RealType pi_times(RealType e, RealType& lo) const {
-        RealType p = pi_hi * e;
-        RealType q = pi_lo * e;
+        RealType p, ep, q, eq;
+        product(pi_hi, e, p, ep);
+        product(pi_lo, e, q, eq);
         RealType hi = p + q;
         // |q| << |p|, so the sum's error is simply q - (hi - p). The rounding
         // of q itself only matters for the large multiples of the argument
         // reduction, but it is cheap to carry.
-        lo = (q - (hi - p)) + jacobi_theta_product_error(pi_hi, e, p, splitter) + jacobi_theta_product_error(pi_lo, e, q, splitter);
+        lo = (q - (hi - p)) + ep + eq;
         return hi;
     }
 
     // (N + dN) / (a_hi + a_lo) as E + dE, to second order
     RealType divide_by_a(RealType N, RealType dN, RealType& dE) const {
         RealType E = N * inv_a_hi;
-        RealType Ea = E * a_hi;
-        RealType r = (N - Ea) - jacobi_theta_product_error(E, a_hi, Ea, splitter);
+        RealType Ea, eEa;
+        product(E, a_hi, Ea, eEa);
+        RealType r = (N - Ea) - eEa;
         dE = (r + dN - E * a_lo) * inv_a_hi;
         return E;
     }
@@ -409,7 +335,7 @@ struct jacobi_theta_exponents {
         RealType result = exp(-E);
         if (result == 0)
             return result;
-        RealType dE = jacobi_theta_product_error(a_hi, e, E, splitter) + a_lo * e;
+        RealType dE = product.error(a_hi, e, E) + a_lo * e;
         return result * (1 - dE);
     }
 
@@ -451,12 +377,12 @@ struct jacobi_theta_exponents {
             RealType c = pi_times(e, dc);
             dc += pi_lo2 * e;
             RealType r = reduced - c; // exact, c being within a factor of two of reduced
-            RealType err = (jacobi_theta_sum_error(reduced, RealType(-c), r) - dc) + dz;
+            RealType err = (two_sum_error(reduced, RealType(-c), r) - dc) + dz;
             // err can be as large as half an ulp of z, too large for the
             // first-order corrections applied to the terms, so it is folded
             // into the remainder, leaving only the rounding of that.
             reduced = r + err;
-            dz = jacobi_theta_sum_error(r, err, reduced);
+            dz = two_sum_error(r, err, reduced);
             // m may exceed the integers the type holds exactly, but such
             // values are even, so their parity is still right.
             k += fmod(m, RealType(2));
@@ -473,10 +399,11 @@ struct jacobi_theta_exponents {
             RealType dc;
             RealType c = pi_times(RealType(m), dc) / 2;
             z_n = z + c;
-            dz += jacobi_theta_sum_error(z, c, z_n) + dc / 2;
+            dz += two_sum_error(z, c, z_n) + dc / 2;
         }
-        RealType s = z_n * z_n;
-        RealType ds = jacobi_theta_product_error(z_n, z_n, s, splitter) + 2 * z_n * dz;
+        RealType s, ds;
+        product(z_n, z_n, s, ds);
+        ds += 2 * z_n * dz;
         RealType dE;
         RealType E = divide_by_a(s, ds, dE);
         RealType result = exp(-E);
@@ -489,8 +416,9 @@ struct jacobi_theta_exponents {
     // k, where dz is the error in z
     RealType expm1_scaled(RealType z, RealType dz, RealType k) const {
         BOOST_MATH_STD_USING
-        RealType P = z * k;
-        RealType dP = jacobi_theta_product_error(z, k, P, splitter) + dz * k;
+        RealType P, dP;
+        product(z, k, P, dP);
+        dP += dz * k;
         RealType dN;
         RealType N = pi_times(P, dN);
         dN += pi_hi * dP;
@@ -503,7 +431,7 @@ struct jacobi_theta_exponents {
     }
 
 private:
-    jacobi_theta_exponents() : splitter(jacobi_theta_splitter<RealType>()) {
+    jacobi_theta_exponents() {
         pi_hi = jacobi_theta_pi(pi_lo, pi_lo2);
     }
 };
diff --git a/include/boost/math/special_functions/laguerre.hpp b/include/boost/math/special_functions/laguerre.hpp
index e349b5d1ed..19db853d41 100644
--- a/include/boost/math/special_functions/laguerre.hpp
+++ b/include/boost/math/special_functions/laguerre.hpp
@@ -14,7 +14,7 @@
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/tools/config.hpp>
 #include <boost/math/policies/error_handling.hpp>
-#include <boost/math/special_functions/detail/two_product.hpp>
+#include <boost/math/tools/compensated_arithmetic.hpp>
 
 BOOST_MATH_NAMESPACE_BEGIN
 
@@ -46,24 +46,21 @@ T laguerre_imp(unsigned n, T x)
 
    // Errors in p0 and p1, starting with that of 1 - x:
    T e0 = 0;
-   T t = p1 - 1;
-   T e1 = (1 - (p1 - t)) + (-x - t);
+   T e1 = two_sum_error(T(1), T(-x), p1);
    detail::two_product<T> product;
 
    for(unsigned c = 1; c < n; ++c)
    {
       // a + ea = 2c + 1 - x exactly:
       T a = (2 * c + 1) - x;
-      t = a - (2 * c + 1);
-      T ea = ((2 * c + 1) - (a - t)) + (-x - t);
+      T ea = two_sum_error(T(2 * c + 1), T(-x), a);
       // s1 + r1 = a p1 and s2 + r2 = c p0 exactly:
       T s1, r1, s2, r2;
       product(a, p1, s1, r1);
       product(T(c), p0, s2, r2);
       // d + ed = s1 - s2 exactly:
       T d = s1 - s2;
-      t = d - s1;
-      T ed = (s1 - (d - t)) + (-s2 - t);
+      T ed = two_sum_error(s1, T(-s2), d);
       T p2 = d / (c + 1);
       // d - (c + 1) p2, where d - q is exact as q is so close to d:
       T q, rq;
diff --git a/include/boost/math/tools/compensated_arithmetic.hpp b/include/boost/math/tools/compensated_arithmetic.hpp
new file mode 100644
index 0000000000..57b566cbca
--- /dev/null
+++ b/include/boost/math/tools/compensated_arithmetic.hpp
@@ -0,0 +1,182 @@
+//  (C) Copyright Nick Thompson 2026.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_TOOLS_COMPENSATED_ARITHMETIC_HPP
+#define BOOST_MATH_TOOLS_COMPENSATED_ARITHMETIC_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/tools/config.hpp>
+#include <boost/math/tools/precision.hpp>
+
+#if !defined(BOOST_MATH_BUILD_MODULE) && !defined(BOOST_MATH_HAS_NVRTC)
+#include <cmath>
+#endif
+
+//
+// Reassociation (-ffast-math, /fp:fast) simplifies every error term below to zero, so forbid it here.  GCC won't
+// inline these into fast-math callers, so only pin it there when it's needed.  Clang has the pragma from 12
+// (Apple clang 13).
+//
+#if defined(__clang__) && (__clang_major__ >= (defined(__apple_build_version__) ? 13 : 12))
+#define BOOST_MATH_COMPENSATED_ARITHMETIC_CLANG_FP
+#endif
+#if defined(BOOST_MATH_COMPENSATED_ARITHMETIC_CLANG_FP)
+#pragma float_control(push)
+#pragma clang fp reassociate(off)
+#elif defined(__GNUC__) && defined(__ASSOCIATIVE_MATH__)
+#pragma GCC push_options
+#pragma GCC optimize("no-associative-math")
+#elif defined(_MSC_VER) && !defined(__clang__) && defined(_M_FP_FAST)
+#pragma float_control(precise, on, push)
+#endif
+
+BOOST_MATH_NAMESPACE_BEGIN namespace detail {
+
+#ifndef BOOST_MATH_HAS_NVRTC
+
+namespace two_product_fma {
+using std::fma;
+
+// result = fma(a, b, c) where the type has an fma, found by argument-dependent lookup for the
+// non-built-in types.  Returns false, leaving result alone, where there is none.
+template <class T>
+inline auto evaluate(const T& a, const T& b, const T& c, T& result, int) -> decltype(static_cast<void>(fma(a, b, c)), bool())
+{
+   result = fma(a, b, c);
+   return true;
+}
+template <class T>
+inline bool evaluate(const T&, const T&, const T&, T&, long)
+{
+   return false;
+}
+} // namespace two_product_fma
+
+//
+// Splits a product into its rounded value and its rounding error: p + e = a b exactly.  A fused fma gives e
+// directly, but having an fma is no guarantee of one: some C libraries (newlib, hence Cygwin) and the
+// Boost.Multiprecision backends without a native one (cpp_bin_float) evaluate it as a * b + c with two
+// roundings, which would make every e zero.  So probe once per type: (1 + eps)^2 = 1 + 2 eps + eps^2, and only
+// a fused fma(x, x, -(1 + 2 eps)) returns the last term.  Otherwise split the factors in half (Veltkamp) so that
+// the partial products are exact (Dekker).
+// The rounded product comes from fma(a, b, 0), which the compiler can't contract into a later addition.
+//
+template <class T>
+class two_product
+{
+public:
+   two_product() : splitter(0)
+   {
+      BOOST_MATH_STD_USING
+      if (!is_fused())
+         splitter = ldexp(T(1), (tools::digits<T>() + 1) / 2) + 1;
+   }
+   void operator()(const T& a, const T& b, T& p, T& e) const
+   {
+      if (is_fused())
+         two_product_fma::evaluate(a, b, T(0), p, 0);
+      else
+         p = a * b;
+      e = error(a, b, p);
+   }
+   // The rounding error of ab = fl(a b): a b = ab + error(a, b, ab) exactly.
+   T error(const T& a, const T& b, const T& ab) const
+   {
+      if (is_fused())
+      {
+         T e;
+         two_product_fma::evaluate(a, b, T(-ab), e, 0);
+         return e;
+      }
+      T ah, al, bh, bl;
+      split(a, ah, al);
+      split(b, bh, bl);
+      return (((ah * bh - ab) + ah * bl) + al * bh) + al * bl;
+   }
+private:
+   void split(const T& t, T& hi, T& lo) const
+   {
+      T g = splitter * t;
+      hi = g - (g - t);
+      lo = t - hi;
+   }
+   static bool is_fused()
+   {
+      static const bool result = []
+      {
+         // From a volatile, so that the compiler can't fold the probe at compile time:
+         volatile int volatile_one = 1;
+         T one = static_cast<int>(volatile_one);
+         T eps = tools::epsilon<T>();
+         T x = one + eps;
+         T eps_squared = 0;
+         return two_product_fma::evaluate(x, x, T(-(one + 2 * eps)), eps_squared, 0) && (eps_squared != 0);
+      }();
+      return result;
+   }
+   T splitter;
+};
+
+//
+// a b - c d without the cancellation: with p1 + e1 = a b and p2 + e2 = c d exactly, p1 - p2 is exact where they
+// cancel (Sterbenz's lemma).  Measured near cancellation it was correctly rounded, where Kahan's
+// fma(a, b, -p2) - e2 is out by up to 1.5 ulps.
+//
+template <class T>
+inline T difference_of_products(const T& a, const T& b, const T& c, const T& d)
+{
+   two_product<T> product;
+   T p1, e1, p2, e2;
+   product(a, b, p1, e1);
+   product(c, d, p2, e2);
+   return (p1 - p2) + (e1 - e2);
+}
+
+#endif // BOOST_MATH_HAS_NVRTC
+
+//
+// The exact sum a + b is s + (return value), where s = fl(a + b):
+//
+template <class T>
+BOOST_MATH_GPU_ENABLED inline T two_sum_error(T a, T b, T s)
+{
+   T bb = s - a;
+   return (a - (s - bb)) + (b - bb);
+}
+
+//
+// Neumaier's compensated sum: sum + comp carries the total to about one rounding, however many terms are added.
+//
+template <class T>
+struct compensated_sum
+{
+   T sum;
+   T comp;
+   BOOST_MATH_GPU_ENABLED explicit compensated_sum(T s0 = 0) : sum(s0), comp(0) {}
+   BOOST_MATH_GPU_ENABLED void add(T x)
+   {
+      BOOST_MATH_STD_USING
+      T t = sum + x;
+      comp += (fabs(sum) >= fabs(x)) ? T((sum - t) + x) : T((x - t) + sum);
+      sum = t;
+   }
+   BOOST_MATH_GPU_ENABLED T value() const { return sum + comp; }
+};
+
+} BOOST_MATH_NAMESPACE_END // namespace detail
+
+#if defined(BOOST_MATH_COMPENSATED_ARITHMETIC_CLANG_FP)
+#undef BOOST_MATH_COMPENSATED_ARITHMETIC_CLANG_FP
+#pragma float_control(pop)
+#elif defined(__GNUC__) && defined(__ASSOCIATIVE_MATH__)
+#pragma GCC pop_options
+#elif defined(_MSC_VER) && !defined(__clang__) && defined(_M_FP_FAST)
+#pragma float_control(pop)
+#endif
+
+#endif // BOOST_MATH_TOOLS_COMPENSATED_ARITHMETIC_HPP