From 2939f8fd69f1d26216230d692dc71782504c57e4 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Wed, 30 Sep 2026 13:44:48 +0300 Subject: [PATCH 1/3] Re-express the Chebyshev weights to avoid cancellation All three Chebyshev rules form a quantity that cancels as the nodes approach the interval ends, and in each case the node is a cosine, so the cancelling quantity is available exactly: GaussChebyshev1 w = (pi/n) sqrt(1-x^2) x = cos t -> (pi/n) sin t GaussChebyshev2 w = h sqrt(1-x^2) x = cos t -> h sin t GaussChebyshev3 w = 2h x sqrt((1-x)/x) x = cos^2 t -> h sin(2t) The third collapses further than the first two: 2h cos^2(t) tan(t) is 2h cos(t) sin(t), so the whole product is h sin(2t). Worst-case relative error in the weights against a 50-digit reference: n=50 n=200 n=800 GaussChebyshev1 5.55e-14 -> 2.15e-15 9.44e-13 -> 4.65e-14 1.78e-11 -> 1.13e-13 GaussChebyshev2 5.87e-15 -> 1.67e-15 9.02e-14 -> 1.05e-14 2.93e-12 -> 1.34e-13 GaussChebyshev3 2.44e-13 -> 6.86e-15 6.74e-13 -> 2.08e-14 2.81e-11 -> 4.33e-14 The error grows as O(n^2) before the change because 1 - cos^2 t is O(1/n^2) against terms of O(1), so it is worst exactly where a large radial grid is most likely to be trusted. Two honest caveats about what this does and does not buy. It changes no integral at double precision. These rules have their weights rescaled by sqrt(1-x^2), so they are not Gauss rules for polynomials and their truncation error dominates: integrating 1 and x^2 gives identical results before and after, at 1e-4 to 1e-7. The gain is the fidelity of the weights to their defining formula, which matters for a consumer that amplifies them -- a radial Jacobian of 1/(1-x) or 1/(1-x)^2 -- and for anything comparing against a reference. It is not a double-precision artefact. The loss is relative and so independent of working precision: at 40 digits, GaussChebyshev1 at n=800 still loses two orders, 5.62e-37 against 4.11e-39 re-expressed. Any arbitrary-precision use of these rules inherits it. The existing tests pass either way and could not have caught this: they check weights with an absolute tolerance of 10*eps, and the weights near the endpoints are O(1e-4), so the effective relative tolerance is ~2e-11. Two further sites from the same audit are left open, both needing a design decision rather than a local edit: the radial transforms' own 1-x, which the transform cannot re-express because only the primitive knows the node is a cosine, and the GaussChebyshev2Modified node, which has no exact re-expression and needs a series. See #118 and #117. Refs #118. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- .../quadratures/primitive/gausschebyshev1.hpp | 7 +++++-- .../quadratures/primitive/gausschebyshev2.hpp | 6 ++++-- .../quadratures/primitive/gausschebyshev3.hpp | 12 ++++++++---- 3 files changed, 17 insertions(+), 8 deletions(-) diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp index f07427e..fc371e1 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp @@ -77,8 +77,11 @@ struct quadrature_traits> { 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). + wi *= std::sin(ti); // Store into memory points[idx] = xi; diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp index ef64afe..46ef4b2 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp @@ -78,8 +78,10 @@ struct quadrature_traits> { 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. + const auto wi = pi_ov_npts_p_1 * std::sin(ti); // Store into memory points[idx] = xi; diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp index bb242d8..dda5135 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp @@ -72,11 +72,15 @@ struct quadrature_traits> { 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. + const auto wi = pi_ov_2n_p_1 * std::sin(2.0 * ti); // Copy to storage points[idx] = xi; From 030db773c61245c086b6f469000a2ca068955ab0 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Tue, 6 Oct 2026 09:35:03 +0300 Subject: [PATCH 2/3] Evaluate the Chebyshev sines at the complementary angle Re-expressing the weights as sin(t) removed the 1 - x*x cancellation but left a second, larger error in place. All three generators walk the index as i = npts - idx, so t lies in (pi/2,pi), where sin(t) ~ pi - t. The O(ulp(pi)) absolute rounding of t is then an O(n eps) relative error in the weight, and the reflection copies it to both endpoints. The complement is exact in the loop index, so evaluate the sine there instead: (2 idx + 1) h in T1, (idx + 1) h in T2, and min(2i-1, 2n+2-2i) h in T3, which reflects 2t about pi/2. Worst relative weight error against a 50-digit reference at n=800, by rule: before after T1 1.1e-13 3.6e-16 T2 1.3e-13 4.7e-16 T3 4.3e-14 2.8e-16 The remaining error is flat in n rather than growing, and the nodes are untouched. The new test case pins the extremal weights at n=800 to a relative 1e-14, which the existing absolute 10*eps checks cannot see, since the weights are themselves O(1/n^2) there. Its reference is the Taylor series of sin about zero, independent of std::sin and accurate to rounding at these angles. Five of its six assertions fail without this change. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- .../quadratures/primitive/gausschebyshev1.hpp | 8 ++- .../quadratures/primitive/gausschebyshev2.hpp | 6 ++- .../quadratures/primitive/gausschebyshev3.hpp | 8 ++- test/1d_quadratures.cxx | 52 +++++++++++++++++++ 4 files changed, 71 insertions(+), 3 deletions(-) diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp index fc371e1..4238886 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp @@ -81,7 +81,13 @@ struct quadrature_traits> { // 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). - wi *= std::sin(ti); + // + // 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; diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp index 46ef4b2..335e6ac 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp @@ -81,7 +81,11 @@ struct quadrature_traits> { // [-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. - const auto wi = pi_ov_npts_p_1 * std::sin(ti); + // + // 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; diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp index dda5135..3d3b35c 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp @@ -80,7 +80,13 @@ struct quadrature_traits> { // 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. - const auto wi = pi_ov_2n_p_1 * std::sin(2.0 * ti); + // + // 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; diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index a7cc128..7e5ce3e 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -184,6 +184,58 @@ 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 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 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") { + // w = h sin(2 h) at the node nearest 1, and h sin(h) at the one nearest 0 + const double h = M_PI / (2.0 * npts + 1); + const GaussChebyshev3 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 quad(150); const auto msg = "Euler-Maclaurin N = " + std::to_string(quad.npts()); From 70c0ffe98e2755f2963420f60d7322df66857cd5 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Fri, 9 Oct 2026 10:15:08 +0300 Subject: [PATCH 3/3] Fix the swapped endpoint labels in the T3 comment In the third-kind rule the loop walks i = npts - idx, so front() is the node nearest 0, where t = pi/2 - h exactly and x = sin^2(h), carrying h sin(2h); back() is the node nearest 1, carrying h sin(h). The comment had the two the wrong way round. The assertions were already correct. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- test/1d_quadratures.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 7e5ce3e..fcc5591 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -226,7 +226,9 @@ TEST_CASE( "Gauss-Chebyshev extremal weights at large N", "[1d-quad]" ) { } SECTION("T3") { - // w = h sin(2 h) at the node nearest 1, and h sin(h) at the one nearest 0 + // 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 q(npts); const auto& wgts = q.weights();