Keep ibeta monotonic under -ffp-contract=fast and fix develop's CI

.github/workflows/ci.yml

@@ -14,8 +14,9 @@ on:
   release:
     types: [published, created, edited]
 
+# A new push to a branch (develop, or a PR's) cancels that branch's older runs of this workflow:
 concurrency:
-  group: ${{ github.head_ref || github.run_id }}
+  group: ${{ github.workflow }}-${{ github.head_ref || github.ref }}
   cancel-in-progress: true
 
 # Every job builds with the same floating-point contraction: fused multiply-add, contracted across

.github/workflows/codecov.yml

@@ -19,6 +19,11 @@ on:
       - fix/**
       - pr/**
 
+# A new push to a branch (develop, or a PR's) cancels that branch's older runs of this workflow:
+concurrency:
+  group: ${{ github.workflow }}-${{ github.head_ref || github.ref }}
+  cancel-in-progress: true
+
 env:
   GIT_FETCH_JOBS: 8
   NET_RETRY_COUNT: 5

.gitignore

@@ -32,3 +32,7 @@ build.ninja
 a.out
 *.json
 
+# Test executables built in place (extensionless), but not the sources or the test_* directories:
+/test/test_*
+!/test/test_*.*
+!/test/test_*/

include/boost/math/differentiation/lanczos_smoothing.hpp

@@ -50,6 +50,9 @@ class discrete_legendre {
         return prod;
     }
 
+    // Where a step adds two products, contraction may fuse either one, and differently in each inlined copy.
+    // fma(a, b, 0) rounds one of them on its own, which leaves the compiler a single choice; types that don't
+    // contract (x87 long double) round exactly as before.
     Real next()
     {
         Real N = 2 * m_n + 1;
@@ -63,10 +66,11 @@ class discrete_legendre {
 
     Real next_prime()
     {
+        using std::fma;
         Real N = 2 * m_n + 1;
         Real s = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n);
-        Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - s * m_qrm2) / m_r;
-        Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - s * m_qrm2p) / m_r;
+        Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - fma(s, m_qrm2, Real(0))) / m_r;
+        Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - fma(s, m_qrm2p, Real(0))) / m_r;
         m_qrm2 = m_qrm1;
         m_qrm1 = tmp1;
         m_qrm2p = m_qrm1p;
@@ -77,12 +81,13 @@ class discrete_legendre {
 
     Real next_dbl_prime()
     {
+        using std::fma;
         Real N = 2*m_n + 1;
         Real trm1 = 2*m_r - 1;
         Real s = (m_r - 1) * (N * N - (m_r - 1) * (m_r - 1)) / Real(4 * m_n * m_n);
-        Real rqrpp = 2*trm1*m_qrm1p + trm1*m_x*m_qrm1pp - s*m_qrm2pp;
-        Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - s * m_qrm2) / m_r;
-        Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - s * m_qrm2p) / m_r;
+        Real rqrpp = fma(2*trm1, m_qrm1p, Real(0)) + trm1*m_x*m_qrm1pp - fma(s, m_qrm2pp, Real(0));
+        Real tmp1 = ((2 * m_r - 1) * m_x * m_qrm1 - fma(s, m_qrm2, Real(0))) / m_r;
+        Real tmp2 = ((2 * m_r - 1) * (m_qrm1 + m_x * m_qrm1p) - fma(s, m_qrm2p, Real(0))) / m_r;
         m_qrm2 = m_qrm1;
         m_qrm1 = tmp1;
         m_qrm2p = m_qrm1p;
@@ -93,53 +98,32 @@ class discrete_legendre {
         return m_qrm1pp;
     }
 
-    Real operator()(Real x, std::size_t k)
+    Real operator()(Real x, std::size_t k) const
     {
         BOOST_MATH_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required.");
         if (k == 0)
         {
             return 1;
         }
-        if (k == 1)
-        {
-            return x;
-        }
-        Real qrm2 = 1;
-        Real qrm1 = x;
-        Real N = 2 * m_n + 1;
+        discrete_legendre dlp(m_n, x);
+        Real q = x;
         for (std::size_t r = 2; r <= k; ++r) {
-            Real num = (r - 1) * (N * N - (r - 1) * (r - 1)) * qrm2;
-            Real tmp = (2 * r - 1) * x * qrm1 - num / Real(4 * m_n * m_n);
-            qrm2 = qrm1;
-            qrm1 = tmp / r;
+            q = dlp.next();
         }
-        return qrm1;
+        return q;
     }
 
-    Real prime(Real x, std::size_t k) {
+    Real prime(Real x, std::size_t k) const {
         BOOST_MATH_ASSERT_MSG(k <= 2 * m_n, "r <= 2n is required.");
         if (k == 0) {
             return 0;
         }
-        if (k == 1) {
-            return 1;
-        }
-        Real qrm2 = 1;
-        Real qrm1 = x;
-        Real qrm2p = 0;
-        Real qrm1p = 1;
-        Real N = 2 * m_n + 1;
+        discrete_legendre dlp(m_n, x);
+        Real qp = 1;
         for (std::size_t r = 2; r <= k; ++r) {
-            Real s =
-                (r - 1) * (N * N - (r - 1) * (r - 1)) / Real(4 * m_n * m_n);
-            Real tmp1 = ((2 * r - 1) * x * qrm1 - s * qrm2) / r;
-            Real tmp2 = ((2 * r - 1) * (qrm1 + x * qrm1p) - s * qrm2p) / r;
-            qrm2 = qrm1;
-            qrm1 = tmp1;
-            qrm2p = qrm1p;
-            qrm1p = tmp2;
+            qp = dlp.next_prime();
         }
-        return qrm1p;
+        return qp;
     }
 
   private:

include/boost/math/interpolators/detail/quintic_hermite_detail.hpp

@@ -140,11 +140,17 @@ class quintic_hermite_detail {
         Real a0 = d2ydx2_[i];
         Real a1 = d2ydx2_[i+1];
         Real t= (x-x0)/dx;
-        Real t2 = t*t;
 
-        Real dydx = 30*t2*(1 - 2*t + t*t)*(y1-y0)/dx;
-        dydx += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*v0 - t*t*(12 - 28*t + 15*t*t)*v1;
-        dydx += (t*dx/2)*((2 - 9*t + 12*t*t - 5*t*t*t)*a0 + t*(3 - 8*t + 5*t*t)*a1);
+        // The basis polynomials factored at their zeros, with u = 1 - t (exact for t >= 1/2), as in septic_hermite:
+        Real u = 1 - t;
+        Real z0 = 30*t*t*u*u;
+        Real z1 = u*u*(1 + t*(2 - 15*t));
+        Real z2 = t*t*(1 + u*(2 - 15*u));
+        // z0 m + z1 v0 + z2 v1 (m the secant slope) cancels; with z0 + z1 + z2 = 1 it is anchored at the
+        // nearer endpoint's derivative instead, as in septic_hermite:
+        Real m = (y1-y0)/dx;
+        Real dydx = t < Real(0.5) ? Real(v0 + (z0*(m - v0) + z2*(v1 - v0))) : Real(v1 + (z0*(m - v1) + z1*(v0 - v1)));
+        dydx += (t*dx/2)*(u*u*(2 - 5*t)*a0 + t*u*(3 - 5*t)*a1);
         return dydx;
     }
 
@@ -334,9 +340,16 @@ class cardinal_quintic_hermite_detail {
         Real d2y0 = d2y_[i];
         Real d2y1 = d2y_[i+1];
 
-        Real dydx = 30*t*t*(1 - 2*t + t*t)*(y1-y0);
-        dydx += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*dy0 - t*t*(12 - 28*t + 15*t*t)*dy1;
-        dydx += t*((2 - 9*t + 12*t*t - 5*t*t*t)*d2y0 + t*(3 - 8*t + 5*t*t)*d2y1);
+        // The basis polynomials factored at their zeros, with u = 1 - t (exact for t >= 1/2), as in septic_hermite:
+        Real u = 1 - t;
+        Real z0 = 30*t*t*u*u;
+        Real z1 = u*u*(1 + t*(2 - 15*t));
+        Real z2 = t*t*(1 + u*(2 - 15*u));
+        // z0 m + z1 v0 + z2 v1 (m the secant slope) cancels; with z0 + z1 + z2 = 1 it is anchored at the
+        // nearer endpoint's derivative instead, as in septic_hermite:
+        Real m = y1 - y0;
+        Real dydx = t < Real(0.5) ? Real(dy0 + (z0*(m - dy0) + z2*(dy1 - dy0))) : Real(dy1 + (z0*(m - dy1) + z1*(dy0 - dy1)));
+        dydx += t*(u*u*(2 - 5*t)*d2y0 + t*u*(3 - 5*t)*d2y1);
         dydx *= inv_dx_;
         return dydx;
     }
@@ -517,9 +530,16 @@ class cardinal_quintic_hermite_detail_aos {
         Real a0 = data_[i][2];
         Real a1 = data_[i+1][2];
 
-        Real dy = 30*t*t*(1 - 2*t + t*t)*(y1-y0);
-        dy += (1-18*t*t + 32*t*t*t - 15*t*t*t*t)*v0 - t*t*(12 - 28*t + 15*t*t)*v1;
-        dy += t*((2 - 9*t + 12*t*t - 5*t*t*t)*a0 + t*(3 - 8*t + 5*t*t)*a1);
+        // The basis polynomials factored at their zeros, with u = 1 - t (exact for t >= 1/2), as in septic_hermite:
+        Real u = 1 - t;
+        Real z0 = 30*t*t*u*u;
+        Real z1 = u*u*(1 + t*(2 - 15*t));
+        Real z2 = t*t*(1 + u*(2 - 15*u));
+        // z0 m + z1 v0 + z2 v1 (m the secant slope) cancels; with z0 + z1 + z2 = 1 it is anchored at the
+        // nearer endpoint's derivative instead, as in septic_hermite:
+        Real m = y1 - y0;
+        Real dy = t < Real(0.5) ? Real(v0 + (z0*(m - v0) + z2*(v1 - v0))) : Real(v1 + (z0*(m - v1) + z1*(v0 - v1)));
+        dy += t*(u*u*(2 - 5*t)*a0 + t*u*(3 - 5*t)*a1);
         return dy*inv_dx_;
     }
 

include/boost/math/interpolators/detail/septic_hermite_detail.hpp

@@ -162,16 +162,24 @@ class septic_hermite_detail {
         Real t = (x-x0)/dx;
         Real t2 = t*t;
         Real t3 = t2*t;
-        Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
-        Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
-        Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
-        Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
-        Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
-        Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
-        Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
-
-        Real dydx = z0*(y1-y0)/dx;
-        dydx += z1*v0 + z2*v1;
+        // The basis polynomials factored at their zeros, with u = 1 - t (exact for t >= 1/2): expanded,
+        // they lose up to ~60 eps to cancellation near their zeros.
+        Real u = 1 - t;
+        Real u2 = u*u;
+        Real u3 = u2*u;
+        Real z1 = u3*(1 + t*(3 + t*(6 - 70*t)));
+        Real z2 = t3*(1 + u*(3 + u*(6 - 70*u)));
+        Real z3 = u3*(1 + t*(3 - 14*t));
+        Real z4 = t2*u*(10 + t*(-25 + 14*t));
+        Real z5 = u3*(3 - 7*t);
+        Real z6 = t*u2*(-4 + 7*t);
+
+        // z0 m + z1 v0 + z2 v1, with m the secant slope, has terms ~10x the result that cancel. Since
+        // z0 + z1 + z2 = 1 (lines are reproduced), it is v0 + z0 (m - v0) + z2 (v1 - v0), whose terms
+        // vanish at t = 0, or, nearer t = 1, the same anchored at v1:
+        Real m = (y1-y0)/dx;
+        Real z0 = 140*t3*u3;
+        Real dydx = t < Real(0.5) ? Real(v0 + (z0*(m - v0) + z2*(v1 - v0))) : Real(v1 + (z0*(m - v1) + z1*(v0 - v1)));
         dydx += (x-x0)*(z3*a0 + z4*a1);
         dydx += (x-x0)*(x-x0)*(z5*j0 + z6*j1)/6;
         return dydx;
@@ -351,16 +359,22 @@ class cardinal_septic_hermite_detail {
         Real j1 = d3y_[i+1];
         Real t2 = t*t;
         Real t3 = t2*t;
-        Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
-        Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
-        Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
-        Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
-        Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
-        Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
-        Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
-
-        Real dydx = z0*(y1-y0)*inv_dx_;
-        dydx += (z1*dy0 + z2*dy1)*inv_dx_;
+        // The basis polynomials factored at their zeros, with u = 1 - t (exact for t >= 1/2): expanded,
+        // they lose up to ~60 eps to cancellation near their zeros.
+        Real u = 1 - t;
+        Real u2 = u*u;
+        Real u3 = u2*u;
+        Real z1 = u3*(1 + t*(3 + t*(6 - 70*t)));
+        Real z2 = t3*(1 + u*(3 + u*(6 - 70*u)));
+        Real z3 = u3*(1 + t*(3 - 14*t));
+        Real z4 = t2*u*(10 + t*(-25 + 14*t));
+        Real z5 = u3*(3 - 7*t);
+        Real z6 = t*u2*(-4 + 7*t);
+
+        // As in septic_hermite_detail::prime, with the slope y1 - y0 in units of dx:
+        Real m = y1 - y0;
+        Real z0 = 140*t3*u3;
+        Real dydx = (t < Real(0.5) ? Real(dy0 + (z0*(m - dy0) + z2*(dy1 - dy0))) : Real(dy1 + (z0*(m - dy1) + z1*(dy0 - dy1))))*inv_dx_;
         dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
         dydx += t*t*(z5*j0 + z6*j1)*inv_dx_;
         return dydx;
@@ -562,16 +576,22 @@ class cardinal_septic_hermite_detail_aos {
         Real j1 = data_[i+1][3];
         Real t2 = t*t;
         Real t3 = t2*t;
-        Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
-        Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
-        Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
-        Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
-        Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
-        Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
-        Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
-
-        Real dydx = z0*(y1-y0)*inv_dx_;
-        dydx += (z1*dy0 + z2*dy1)*inv_dx_;
+        // The basis polynomials factored at their zeros, with u = 1 - t (exact for t >= 1/2): expanded,
+        // they lose up to ~60 eps to cancellation near their zeros.
+        Real u = 1 - t;
+        Real u2 = u*u;
+        Real u3 = u2*u;
+        Real z1 = u3*(1 + t*(3 + t*(6 - 70*t)));
+        Real z2 = t3*(1 + u*(3 + u*(6 - 70*u)));
+        Real z3 = u3*(1 + t*(3 - 14*t));
+        Real z4 = t2*u*(10 + t*(-25 + 14*t));
+        Real z5 = u3*(3 - 7*t);
+        Real z6 = t*u2*(-4 + 7*t);
+
+        // As in septic_hermite_detail::prime, with the slope y1 - y0 in units of dx:
+        Real m = y1 - y0;
+        Real z0 = 140*t3*u3;
+        Real dydx = (t < Real(0.5) ? Real(dy0 + (z0*(m - dy0) + z2*(dy1 - dy0))) : Real(dy1 + (z0*(m - dy1) + z1*(dy0 - dy1))))*inv_dx_;
         dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
         dydx += t*t*(z5*j0 + z6*j1)*inv_dx_;
         return dydx;

include/boost/math/quadrature/gauss_hermite.hpp

@@ -182,6 +182,27 @@ class gauss_hermite_detail
          w[positive_zeros + offset - t] = weight(z, r);
          upper = lower;
       }
+      // Every node runs through the same rounded coefficients a and b, so the weights share a relative error
+      // of a few eps that grows with N (against 50 digits: -5 eps at N = 200, -10 at 500); exact coefficients
+      // remove it only with the recurrence in twice the precision. The exact weights sum to sqrt(pi), the
+      // integral of exp(-x^2), so scale them to that: weights with relative errors e_i then have e_i - e, e
+      // 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;
+      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;
+      }
+      const Real scale = BOOST_MATH_NAMESPACE::constants::root_pi<Real>() / (sum + compensation);
+      for (Real& weight_i : w)
+      {
+         weight_i *= scale;
+      }
       return std::make_pair(x, w);
    }