Declare bound_inclusive on every primitive quadrature

include/integratorxx/quadratures/primitive/gausschebyshev1.hpp

@@ -55,6 +55,8 @@ struct quadrature_traits<GaussChebyshev1<PointType, WeightType>> {
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
 
+  inline static constexpr bool bound_inclusive = false;
+
   inline static std::tuple<point_container, weight_container> generate(size_t npts) {
     point_container points(npts);
     weight_container weights(npts);

include/integratorxx/quadratures/primitive/gausschebyshev2modified.hpp

@@ -59,6 +59,8 @@ struct quadrature_traits<GaussChebyshev2Modified<PointType, WeightType>> {
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
 
+  inline static constexpr bool bound_inclusive = false;
+
   inline static std::tuple<point_container, weight_container> generate(size_t npts) {
     const weight_type oonpp = 1.0 / (npts + 1);
 

include/integratorxx/quadratures/primitive/gausschebyshev3.hpp

@@ -55,6 +55,8 @@ struct quadrature_traits<GaussChebyshev3<PointType, WeightType>> {
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
 
+  inline static constexpr bool bound_inclusive = false;
+
   inline static std::tuple<point_container, weight_container> generate(size_t npts) {
     const weight_type pi_ov_2n_p_1 = M_PI / (2 * npts + 1);
 

include/integratorxx/quadratures/primitive/gausslegendre.hpp

@@ -44,6 +44,8 @@ struct quadrature_traits<
   using point_container  = std::vector< point_type >;
   using weight_container = std::vector< weight_type >;
 
+  inline static constexpr bool bound_inclusive = false;
+
   inline static std::tuple<point_container,weight_container>
   generate( size_t npts ) {
     point_container  points( npts );

include/integratorxx/quadratures/primitive/gausslobatto.hpp

@@ -34,6 +34,8 @@ struct quadrature_traits<GaussLobatto<PointType, WeightType>> {
   using point_container = std::vector<point_type>;
   using weight_container = std::vector<weight_type>;
 
+  inline static constexpr bool bound_inclusive = true;
+
   inline static std::tuple<point_container, weight_container> generate(
       size_t npts) {
     point_container points(npts);

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 <typename Primitive>
+void compose_with_all_transforms() {
+  constexpr size_t npts = 16;
+  CHECK( IntegratorXX::RadialTransformQuadrature<Primitive,
+           IntegratorXX::BeckeRadialTraits>(npts).npts() == npts );
+  CHECK( IntegratorXX::RadialTransformQuadrature<Primitive,
+           IntegratorXX::MuraKnowlesRadialTraits>(npts).npts() == npts );
+  CHECK( IntegratorXX::RadialTransformQuadrature<Primitive,
+           IntegratorXX::MurrayHandyLamingRadialTraits<2>>(npts).npts() == npts );
+  CHECK( IntegratorXX::RadialTransformQuadrature<Primitive,
+           IntegratorXX::TreutlerAhlrichsRadialTraits>(npts).npts() == npts );
+}
+
 TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) {
 
   // Reference integral for polynomial evaluated over [-1,1]
@@ -191,6 +208,50 @@ TEST_CASE( "Becke Quadratures", "[1d-quad]" ) {
   test_quadrature<RadialGaussian>(msg, quad, std::sqrt(M_PI)/4, 1e-10);
 }
 
+
+TEST_CASE( "Radial transforms over the primitive rules", "[1d-quad]" ) {
+
+  // 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 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<double,double>,
+      IntegratorXX::BeckeRadialTraits> quad(350);
+    REQUIRE( quad.npts() == 350 );
+    const auto msg = "Gauss-Lobatto x Becke N = " + std::to_string(quad.npts());
+    test_quadrature<RadialGaussian>(msg, quad, std::sqrt(M_PI)/4, 1e-10);
+  }
+
+  SECTION("Gauss-Legendre x Becke") {
+    IntegratorXX::RadialTransformQuadrature<
+      IntegratorXX::GaussLegendre<double,double>,
+      IntegratorXX::BeckeRadialTraits> quad(350);
+    REQUIRE( quad.npts() == 350 );
+    const auto msg = "Gauss-Legendre x Becke N = " + std::to_string(quad.npts());
+    test_quadrature<RadialGaussian>(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<IntegratorXX::GaussChebyshev1<double,double>>();
+    compose_with_all_transforms<IntegratorXX::GaussChebyshev2<double,double>>();
+    compose_with_all_transforms<IntegratorXX::GaussChebyshev2Modified<double,double>>();
+    compose_with_all_transforms<IntegratorXX::GaussChebyshev3<double,double>>();
+    compose_with_all_transforms<IntegratorXX::GaussLegendre<double,double>>();
+    compose_with_all_transforms<IntegratorXX::GaussLobatto<double,double>>();
+    compose_with_all_transforms<IntegratorXX::UniformTrapezoid<double,double>>();
+  }
+
+}
+
 TEST_CASE( "Lebedev-Laikov", "[1d-quad]" ) {