Re-express the Chebyshev weights to avoid cancellation

include/integratorxx/quadratures/primitive/gausschebyshev1.hpp

@@ -77,8 +77,17 @@ struct quadrature_traits<GaussChebyshev1<PointType, WeightType>> {
       auto wi = pi_ov_npts;
 
       // However, as we're integrating f(x) not \frac{f(x)}{\sqrt{1 -
-      // x^2}}, we must factor the \sqrt{1-x^2} into the weight
-      wi *= std::sqrt(1.0 - xi * xi);
+      // x^2}}, we must factor the \sqrt{1-x^2} into the weight.
+      // Since x_i = cos(t_i) with t_i in (0,pi), sqrt(1-x^2) is sin(t_i)
+      // exactly; forming 1 - x*x instead cancels badly as the nodes approach
+      // +-1, losing accuracy as O(n^2).
+      //
+      // t_i lies in (pi/2,pi) here, where sin(t) ~ pi - t, so the O(ulp(pi))
+      // absolute rounding of t_i would show up as an O(n eps) *relative*
+      // error in sin(t_i), and the reflection below copies it to both ends.
+      // The complement is exact in the index, t_i = pi - (2 idx + 1) h, so
+      // evaluate the sine at that small angle instead.
+      wi *= std::sin((2.0 * idx + 1.0) * pi_ov_2npts);
 
       // Store into memory
       points[idx]  = xi;

include/integratorxx/quadratures/primitive/gausschebyshev2.hpp

@@ -78,8 +78,14 @@ struct quadrature_traits<GaussChebyshev2<PointType, WeightType>> {
       const auto xi = std::cos(ti);
 
       // The quadrature weight, transformed to unit weight factor in
-      // [-1,1] is given by (see comments above for explanation)
-      const auto wi = pi_ov_npts_p_1 * std::sqrt(1.0 - xi * xi);
+      // [-1,1] is given by (see comments above for explanation). Since
+      // x_i = cos(t_i) with t_i in (0,pi), sqrt(1-x^2) is sin(t_i) exactly;
+      // forming 1 - x*x instead cancels badly as the nodes approach +-1.
+      //
+      // As in the T1 rule, t_i lies in (pi/2,pi), where the absolute rounding
+      // of t_i is an O(n eps) relative error in sin(t_i). Use the complement,
+      // which is exact in the index: t_i = pi - (idx + 1) h.
+      const auto wi = pi_ov_npts_p_1 * std::sin((idx + 1.0) * pi_ov_npts_p_1);
 
       // Store into memory
       points[idx] = xi;

include/integratorxx/quadratures/primitive/gausschebyshev3.hpp

@@ -72,11 +72,21 @@ struct quadrature_traits<GaussChebyshev3<PointType, WeightType>> {
       const auto ti = 0.5 * (2 * i - 1) * pi_ov_2n_p_1;
       const auto cti = std::cos(ti);
       const auto xi = cti * cti;  // cos^2(t)
-      auto wi = 2.0 * pi_ov_2n_p_1 * xi;
 
-      // However, since we want the rule with a unit weight factor, we
-      // divide the weights by sqrt(x/(1-x)).
-      wi *= std::sqrt((1.0 - xi) / xi);
+      // The standard weight is 2h x_i, and since we want the rule with a unit
+      // weight factor we divide by sqrt(x/(1-x)). With x_i = cos^2(t_i) and
+      // t_i in (0,pi/2) the whole product collapses exactly:
+      //
+      //   2h cos^2(t) sqrt((1-cos^2 t)/cos^2 t) = 2h cos(t) sin(t) = h sin(2t)
+      //
+      // Forming 1 - x instead cancels badly as the nodes approach 1.
+      //
+      // 2 t_i is m h with m = 2i - 1 an exact integer, and for m h past pi/2
+      // the sine of the supplement is the relatively accurate one, since
+      // sin(s) ~ pi - s there. Reflect about pi/2 using the index alone.
+      const size_t m = 2 * i - 1;
+      const size_t k = (m <= npts) ? m : (2 * npts + 1 - m);
+      const auto wi = pi_ov_2n_p_1 * std::sin(k * pi_ov_2n_p_1);
 
       // Copy to storage
       points[idx]  = xi;

test/1d_quadratures.cxx

@@ -184,6 +184,60 @@ TEST_CASE( "Gauss-Chebyshev T3 Quadratures", "[1d-quad]" ) {
 }
 
 
+TEST_CASE( "Gauss-Chebyshev extremal weights at large N", "[1d-quad]" ) {
+
+  // The extremal weight of each Chebyshev rule is sin(t) at an angle that the
+  // generators reach as pi - t. Rounding t to the nearest double is an
+  // O(ulp(pi)) absolute error there, which is an O(N eps) *relative* error in
+  // the weight: ~1e-13 at N=800, growing with N. Evaluating the sine at the
+  // complement, which is exact in the loop index, keeps it at O(eps).
+  //
+  // The reference is the Taylor series of sin about 0. At these angles
+  // (h < 2e-3) the first omitted term is O(h^7) ~ 1e-22 relative, so the
+  // series is good to rounding and is independent of std::sin.
+  auto sin_series = [](double h) {
+    const double h2 = h * h;
+    return h * (1.0 - h2 / 6.0 * (1.0 - h2 / 20.0 * (1.0 - h2 / 42.0)));
+  };
+
+  // An absolute 10*eps check cannot see this: the weights are themselves
+  // O(1/N^2) here, so the error hides below the tolerance.
+  const double rel_tol = 1e-14;
+  const size_t npts = 800;
+
+  SECTION("T1") {
+    // w_0 = (pi/N) sin(h), h = pi/2N, at the node nearest -1
+    const double h = M_PI / (2.0 * npts);
+    const double ref = (M_PI / npts) * sin_series(h);
+    const GaussChebyshev1<double,double> q(npts);
+    const auto& wgts = q.weights();
+    CHECK_THAT( wgts.front(), Catch::Matchers::WithinRel(ref, rel_tol) );
+    CHECK_THAT( wgts.back(),  Catch::Matchers::WithinRel(ref, rel_tol) );
+  }
+
+  SECTION("T2") {
+    // w_0 = h sin(h), h = pi/(N+1)
+    const double h = M_PI / (npts + 1);
+    const double ref = h * sin_series(h);
+    const GaussChebyshev2<double,double> q(npts);
+    const auto& wgts = q.weights();
+    CHECK_THAT( wgts.front(), Catch::Matchers::WithinRel(ref, rel_tol) );
+    CHECK_THAT( wgts.back(),  Catch::Matchers::WithinRel(ref, rel_tol) );
+  }
+
+  SECTION("T3") {
+    // front() is the node nearest 0: there t = pi/2 - h exactly, so
+    // x = sin^2(h) and the weight is h sin(2h). back() is the node nearest
+    // 1, with weight h sin(h).
+    const double h = M_PI / (2.0 * npts + 1);
+    const GaussChebyshev3<double,double> q(npts);
+    const auto& wgts = q.weights();
+    CHECK_THAT( wgts.front(), Catch::Matchers::WithinRel(h * sin_series(2.0*h), rel_tol) );
+    CHECK_THAT( wgts.back(),  Catch::Matchers::WithinRel(h * sin_series(h),     rel_tol) );
+  }
+
+}
+
 TEST_CASE( "Euler-Maclaurin Quadratures by Murray, Handy, and Laming", "[1d-quad]" ) {
   IntegratorXX::MurrayHandyLaming<double,double> quad(150);
   const auto msg = "Euler-Maclaurin N = " + std::to_string(quad.npts());