From 4fdf32c659b60b8584d8e5bc5727724b84c62678 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Wed, 30 Sep 2026 11:34:27 +0300 Subject: [PATCH 1/3] Declare bound_inclusive on every primitive quadrature RadialTransformQuadrature reads base_quad_traits::bound_inclusive to decide whether to generate two extra nodes and drop the endpoints, but only gausschebyshev2 and uniform declared it. The other five failed to compile when used as a radial base: error: 'bound_inclusive' is not a member of 'base_quad_traits' So a radial transform could be composed only with those two, although the README describes radial quadratures as convolutions of a primitive with a transformation, which reads as a composition matrix. The value is whether the rule places nodes on the interval endpoints: false for the open rules -- gausschebyshev1, gausschebyshev2modified, gausschebyshev3 and gausslegendre -- and true for gausslobatto, which includes both endpoints by construction. gausslobatto is the one that matters beyond bookkeeping: it is the only primitive here whose value is true, and the natural base for a Gauss-Lobatto x uniform-phi product sphere rule. All 28 primitive-transform pairs now compile. The domain-matched ones integrate r^2 exp(-r) over (0, inf) correctly, for instance at 120 points: GaussChebyshev2Modified x Becke 0.00e+00 GaussChebyshev2Modified x TreutlerAhlrichs 6.66e-16 GaussLegendre x Becke 0.00e+00 GaussLobatto x Becke 0.00e+00 GaussChebyshev3 x MurrayHandyLaming 8.88e-16 The mismatched pairs compile but are meaningless, since Becke and Treutler-Ahlrichs map from (-1,1) while Mura-Knowles and Murray-Handy-Laming map from (0,1); the traits member says nothing about domain, only about endpoints. Behaviour-neutral for everything shipped: the four radial quadratures are bit-identical to master over 392 node and weight values, since the two primitives they use already declared the member. Fixes #116. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- .../integratorxx/quadratures/primitive/gausschebyshev1.hpp | 4 ++++ .../quadratures/primitive/gausschebyshev2modified.hpp | 4 ++++ .../integratorxx/quadratures/primitive/gausschebyshev3.hpp | 4 ++++ include/integratorxx/quadratures/primitive/gausslegendre.hpp | 4 ++++ include/integratorxx/quadratures/primitive/gausslobatto.hpp | 4 ++++ 5 files changed, 20 insertions(+) diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp index 9b538aa..ab77e4a 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp @@ -55,6 +55,10 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + /// Whether the rule places nodes on the interval endpoints. + /// Open rule: nodes are cos((2i-1)*pi/2n), which never reaches +-1. + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate(size_t npts) { point_container points(npts); weight_container weights(npts); diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp index 83c8c17..4c169dd 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp @@ -59,6 +59,10 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + /// Whether the rule places nodes on the interval endpoints. + /// Open rule: the nodes are interior for every n. + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate(size_t npts) { const weight_type oonpp = 1.0 / (npts + 1); diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp index 39a8034..2069574 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp @@ -55,6 +55,10 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + /// Whether the rule places nodes on the interval endpoints. + /// Open rule on (0,1): nodes are cos^2(...), never 0 or 1. + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate(size_t npts) { const weight_type pi_ov_2n_p_1 = M_PI / (2 * npts + 1); diff --git a/include/integratorxx/quadratures/primitive/gausslegendre.hpp b/include/integratorxx/quadratures/primitive/gausslegendre.hpp index 755c305..088c1d1 100644 --- a/include/integratorxx/quadratures/primitive/gausslegendre.hpp +++ b/include/integratorxx/quadratures/primitive/gausslegendre.hpp @@ -44,6 +44,10 @@ struct quadrature_traits< using point_container = std::vector< point_type >; using weight_container = std::vector< weight_type >; + /// Whether the rule places nodes on the interval endpoints. + /// Open rule: the Legendre roots are strictly inside (-1,1). + inline static constexpr bool bound_inclusive = false; + inline static std::tuple generate( size_t npts ) { point_container points( npts ); diff --git a/include/integratorxx/quadratures/primitive/gausslobatto.hpp b/include/integratorxx/quadratures/primitive/gausslobatto.hpp index 830e83d..ff96f04 100644 --- a/include/integratorxx/quadratures/primitive/gausslobatto.hpp +++ b/include/integratorxx/quadratures/primitive/gausslobatto.hpp @@ -34,6 +34,10 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; + /// Whether the rule places nodes on the interval endpoints. + /// Lobatto includes both endpoints by construction. + inline static constexpr bool bound_inclusive = true; + inline static std::tuple generate( size_t npts) { point_container points(npts); From 00061bc112a6c9c5f4771c9c3c3cee2c2062f3d4 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Tue, 6 Oct 2026 09:47:05 +0300 Subject: [PATCH 2/3] Cover the primitive-transform pairings, and match the existing style The declarations now read like the two that were already in the tree (gausschebyshev2.hpp, uniform.hpp): bare one-liners, with the reason for each value left to generate() just below. The check from the PR description was not in the suite. Nothing instantiated RadialTransformQuadrature over any of these primitives, since the named radial grids are all fixed aliases over GaussChebyshev2, so the capability this unblocks was untested. Add the full 28-pair sweep, which fails to compile without the declarations ("bound_inclusive is not a member of base_quad_traits"), plus r^2 exp(-r) on Gauss-Lobatto x Becke and Gauss-Legendre x Becke. Lobatto is the only bound_inclusive rule among them, so it also exercises the endpoint-dropping path. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- .../quadratures/primitive/gausschebyshev1.hpp | 2 - .../primitive/gausschebyshev2modified.hpp | 2 - .../quadratures/primitive/gausschebyshev3.hpp | 2 - .../quadratures/primitive/gausslegendre.hpp | 2 - .../quadratures/primitive/gausslobatto.hpp | 2 - test/1d_quadratures.cxx | 58 +++++++++++++++++++ 6 files changed, 58 insertions(+), 10 deletions(-) diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp index ab77e4a..f07427e 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp @@ -55,8 +55,6 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - /// Whether the rule places nodes on the interval endpoints. - /// Open rule: nodes are cos((2i-1)*pi/2n), which never reaches +-1. inline static constexpr bool bound_inclusive = false; inline static std::tuple generate(size_t npts) { diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp index 4c169dd..338c6c3 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp @@ -59,8 +59,6 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - /// Whether the rule places nodes on the interval endpoints. - /// Open rule: the nodes are interior for every n. inline static constexpr bool bound_inclusive = false; inline static std::tuple generate(size_t npts) { diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp index 2069574..bb242d8 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp @@ -55,8 +55,6 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - /// Whether the rule places nodes on the interval endpoints. - /// Open rule on (0,1): nodes are cos^2(...), never 0 or 1. inline static constexpr bool bound_inclusive = false; inline static std::tuple generate(size_t npts) { diff --git a/include/integratorxx/quadratures/primitive/gausslegendre.hpp b/include/integratorxx/quadratures/primitive/gausslegendre.hpp index 088c1d1..183fc9b 100644 --- a/include/integratorxx/quadratures/primitive/gausslegendre.hpp +++ b/include/integratorxx/quadratures/primitive/gausslegendre.hpp @@ -44,8 +44,6 @@ struct quadrature_traits< using point_container = std::vector< point_type >; using weight_container = std::vector< weight_type >; - /// Whether the rule places nodes on the interval endpoints. - /// Open rule: the Legendre roots are strictly inside (-1,1). inline static constexpr bool bound_inclusive = false; inline static std::tuple diff --git a/include/integratorxx/quadratures/primitive/gausslobatto.hpp b/include/integratorxx/quadratures/primitive/gausslobatto.hpp index ff96f04..f662248 100644 --- a/include/integratorxx/quadratures/primitive/gausslobatto.hpp +++ b/include/integratorxx/quadratures/primitive/gausslobatto.hpp @@ -34,8 +34,6 @@ struct quadrature_traits> { using point_container = std::vector; using weight_container = std::vector; - /// Whether the rule places nodes on the interval endpoints. - /// Lobatto includes both endpoints by construction. inline static constexpr bool bound_inclusive = true; inline static std::tuple generate( diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index d106801..1ea0121 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -50,6 +50,23 @@ auto chebyshev_T(int n, T x) { + +// Compose one primitive rule with every radial transform. The pairing is only +// legal if the primitive's traits declare bound_inclusive, so instantiating +// these is itself the regression guard for #116. +template +void compose_with_all_transforms() { + constexpr size_t npts = 16; + CHECK( IntegratorXX::RadialTransformQuadrature(npts).npts() == npts ); + CHECK( IntegratorXX::RadialTransformQuadrature(npts).npts() == npts ); + CHECK( IntegratorXX::RadialTransformQuadrature>(npts).npts() == npts ); + CHECK( IntegratorXX::RadialTransformQuadrature(npts).npts() == npts ); +} + TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { // Reference integral for polynomial evaluated over [-1,1] @@ -191,6 +208,47 @@ TEST_CASE( "Becke Quadratures", "[1d-quad]" ) { test_quadrature(msg, quad, std::sqrt(M_PI)/4, 1e-10); } + +TEST_CASE( "Radial transforms over the primitive rules", "[1d-quad]" ) { + + // The named radial grids are all fixed aliases over GaussChebyshev2, so + // nothing above pairs a transform with any other primitive by hand. + + SECTION("Gauss-Lobatto x Becke") { + // Lobatto is the only bound_inclusive primitive, so this is also the + // endpoint-dropping path: the base rule is asked for npts+2 nodes and the + // two that map to r = 0 and r = infinity are skipped. + IntegratorXX::RadialTransformQuadrature< + IntegratorXX::GaussLobatto, + IntegratorXX::BeckeRadialTraits> quad(350); + REQUIRE( quad.npts() == 350 ); + const auto msg = "Gauss-Lobatto x Becke N = " + std::to_string(quad.npts()); + test_quadrature(msg, quad, std::sqrt(M_PI)/4, 1e-10); + } + + SECTION("Gauss-Legendre x Becke") { + IntegratorXX::RadialTransformQuadrature< + IntegratorXX::GaussLegendre, + IntegratorXX::BeckeRadialTraits> quad(350); + REQUIRE( quad.npts() == 350 ); + const auto msg = "Gauss-Legendre x Becke N = " + std::to_string(quad.npts()); + test_quadrature(msg, quad, std::sqrt(M_PI)/4, 1e-10); + } + + SECTION("Every primitive composes with every transform") { + // All 28 pairings. Accuracy is not the point here: several of these are + // not sensible rules, but all of them must build. + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + compose_with_all_transforms>(); + } + +} + TEST_CASE( "Lebedev-Laikov", "[1d-quad]" ) { From ef50142ebc8963f5292edb3f6244512f88a0454d Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Fri, 9 Oct 2026 10:14:59 +0300 Subject: [PATCH 3/3] Correct what the new test case claims about the aliases Becke and Treutler-Ahlrichs are aliases over GaussChebyshev2, but Mura-Knowles and MHL are over UniformTrapezoid, which is bound_inclusive itself, so the existing MK and MHL cases already exercise endpoint dropping. Lobatto is the only bound_inclusive Gauss rule, not the only bound_inclusive primitive. Comments only; the cases are unchanged. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF --- test/1d_quadratures.cxx | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 1ea0121..a7cc128 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -211,13 +211,16 @@ TEST_CASE( "Becke Quadratures", "[1d-quad]" ) { TEST_CASE( "Radial transforms over the primitive rules", "[1d-quad]" ) { - // The named radial grids are all fixed aliases over GaussChebyshev2, so - // nothing above pairs a transform with any other primitive by hand. + // Becke and Treutler-Ahlrichs are fixed aliases over GaussChebyshev2, and + // Mura-Knowles and MHL over UniformTrapezoid, so nothing above pairs a + // transform with any other primitive by hand. SECTION("Gauss-Lobatto x Becke") { - // Lobatto is the only bound_inclusive primitive, so this is also the - // endpoint-dropping path: the base rule is asked for npts+2 nodes and the - // two that map to r = 0 and r = infinity are skipped. + // Lobatto is the only bound_inclusive Gauss rule, so this is the + // endpoint-dropping path for those: the base rule is asked for npts+2 + // nodes and the two that map to r = 0 and r = infinity are skipped. + // UniformTrapezoid is bound_inclusive as well, so the Mura-Knowles and + // MHL cases above already cover that path. IntegratorXX::RadialTransformQuadrature< IntegratorXX::GaussLobatto, IntegratorXX::BeckeRadialTraits> quad(350);