Consolidate the compensated arithmetic into one header

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());

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)
    {

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)

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;
 

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());

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()));

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;

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)
 
 //