Partially Addresses #100

.github/workflows/build_and_test_compiler_zoo.yml

@@ -95,6 +95,44 @@ jobs:
       shell: bash
       run: cmake --build ${{runner.workspace}}/build --target test
 
+  string_reals_build:
+    name: String Reals Build and Test
+    runs-on: ubuntu-latest
+    container:
+        image: dbwy/chemistry
+    strategy:
+      matrix:
+        header_only: [ON,OFF]
+
+    steps:
+    - uses: actions/checkout@v3
+
+    - name: Setup Compiler
+      shell: bash
+      run: $GITHUB_WORKSPACE/.github/workflows/scripts/compiler_setup.sh gnu 12
+
+    - name: Setup Build Type
+      shell: bash
+      run: echo "set(CMAKE_BUILD_TYPE Release CACHE BOOL \"\" FORCE)" >>
+             ${GITHUB_WORKSPACE}/${GH_ACTIONS_TOOLCHAIN}
+
+    - name: Configure CMake
+      shell: bash
+      run: cmake -S $GITHUB_WORKSPACE -B ${{runner.workspace}}/build
+                 -DINTEGRATORXX_ENABLE_STRING_REALS=ON
+                 -DINTEGRATORXX_HEADER_ONLY=${{matrix.header_only}}
+                 -DCMAKE_INSTALL_PREFIX=${{runner.workspace}}/install
+                 -DCMAKE_PREFIX_PATH=${ENV_PREFIX_PATH}
+                 -DCMAKE_TOOLCHAIN_FILE=${GITHUB_WORKSPACE}/${GH_ACTIONS_TOOLCHAIN}
+
+    - name: Build
+      shell: bash
+      run: cmake --build ${{runner.workspace}}/build -j2
+
+    - name: Test
+      shell: bash
+      run: cmake --build ${{runner.workspace}}/build --target test
+
   subproject_build:
     name: Build as Subproject
     needs: release_build

.gitignore

@@ -1 +1,2 @@
 *.*.swp
+build*/

CMakeLists.txt

@@ -4,6 +4,8 @@ cmake_minimum_required( VERSION 3.17 ) # Require CMake 3.17+
 project( IntegratorXX VERSION 1.2.0 LANGUAGES CXX )
 
 option(INTEGRATORXX_HEADER_ONLY "Force header-only build" OFF)
+option(INTEGRATORXX_ENABLE_STRING_REALS
+       "Carry non-integral literals as their decimal source text" OFF)
 
 if(INTEGRATORXX_HEADER_ONLY)
   add_library( integratorxx INTERFACE )
@@ -21,6 +23,14 @@ target_include_directories( integratorxx
     $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>
 )
 
+# Non-integral literals carried as text rather than as a pre-rounded double.
+# Applied to the target rather than to a directory so that, in the non
+# header-only build, the library's own TUs in src/ and its consumers agree on
+# what `ixx_real` is; they would otherwise disagree across the ABI.
+if(INTEGRATORXX_ENABLE_STRING_REALS)
+  target_compile_definitions( integratorxx ${INTEGRATORXX_TARGET_TYPE} ENABLE_STRING_REALS )
+endif()
+
 include(CheckCXXCompilerFlag)
 check_cxx_compiler_flag("-Wno-missing-braces" INTEGRATORXX_HAS_NO_MISSING_BRACES )
 if( INTEGRATORXX_HAS_NO_MISSING_BRACES )

README.md

@@ -111,7 +111,28 @@ the following radial pruning schemes:
 | Robust    | The Psi4 "robust" pruning scheme     | `PruningScheme::Robust`   |
 | Treutler  | The Treutler-Ahlrichs pruning scheme | `PruntinScheme::Treutler` |
 
-
+### Molecular Grids
+
+`include/integratorxx/molecular_grid/` assembles a multi-atom integration
+grid from per-atom atomic grids (as produced by `SphericalGridFactory`) and
+atomic positions, via `IntegratorXX::MolecularGrid`. Cheap, non-materializing
+metadata (point counts, batch counts, bounding boxes) is available
+immediately after construction; the actual point/weight data is only copied
+for batches a caller explicitly asks for, either by atom subset
+(`batches_for_atoms`) or by a global point-index range
+(`batches_for_point_range`), or for the whole molecule at once via the
+lazily-cached `points()`/`weights()`.
+
+`include/integratorxx/molecular_grid/defaults.hpp` supplies element-indexed
+defaults (`MolecularGridDefaults`, `make_atom_instances`) -- per-element
+radius/scaling-factor tables and Fine/UltraFine/SuperFine/GM3/GM5 size
+presets -- for callers that don't want to hand-pick radial/angular sizes per
+element.
+
+`include/integratorxx/molecular_grid/partition_weights.hpp` provides a
+serial-CPU reference implementation of Becke, Stratmann-Scuseria-Frisch
+(SSF), and Laqua-Kussmann-Ochsenfeld (LKO) fuzzy-cell partition weighting via
+`MolecularGrid::apply_partition_weights`.
 
 ## Example Usage
 
@@ -162,6 +183,26 @@ To use the runtime generator header-only, one needs to include
 `<integratorxx/generators/impl/impl.hpp>` **exactly once** per project,
 otherwise duplicate / incompatible symbols will occur.
 
+### High-precision literals
+
+The tabulated solid-angle grids (Lebedev-Laikov, Delley, Ahrens-Beylkin,
+Womersley) are stored as `IntegratorXX::ixx_real`, which is `double` by default.
+A `double` literal has already been rounded by the time any code can inspect it,
+so a type more precise than `double` -- or one that must bound its own error --
+cannot recover the value the table intended.
+
+Setting `INTEGRATORXX_ENABLE_STRING_REALS=ON` makes `ixx_real` a
+`std::string_view` instead, so the tables carry the exact decimal source text of
+each entry and the conversion to the quadrature's value type happens on read.
+The default `IntegratorXX::fp_traits::from_real` still parses that text via
+`double`, so this option changes nothing on its own; it exists so that a type
+which specializes `from_real` can parse the decimal directly. Unless you have
+such a type, leave it `OFF`.
+
+**N.B.** the setting changes the type of a table entry, so it is not ABI-neutral:
+IntegratorXX and everything that includes its headers must be built with the same
+value.
+
 ## Contributing and Bug Reports
 
 We welcome any and all contributions and encourage bug reports. Please use the

include/integratorxx/batch/spherical_micro_batcher.hpp

@@ -2,6 +2,7 @@
 
 #include <integratorxx/composite_quadratures/spherical_quadrature.hpp>
 #include <integratorxx/type_traits.hpp>
+#include <integratorxx/util/fp_traits.hpp>
 #include <chrono>
 #include <iostream>
 #include <iomanip>
@@ -154,13 +155,23 @@ auto partition_box(
 
   std::vector<T> x_part(npart+1), y_part(npart+1), z_part(npart+1);
 
-  const auto delta_x = extent_x / npart;
-  const auto delta_y = extent_y / npart;
-  const auto delta_z = extent_z / npart;
+  // npart/i are raw loop-counter integers, not T -- routed through
+  // fp_traits<T>::from_integer so this compiles for any T (a raw
+  // `size_t * T`/`T / size_t` requires an exact-type overload of
+  // operator*/operator/ that most non-builtin T don't provide, since it's
+  // resolved by template argument deduction, which -- unlike ordinary
+  // overload resolution -- does not apply an implicit conversion to find
+  // one).
+  using traits = fp_traits<T>;
+  const auto npart_T = traits::from_integer(static_cast<ixx_int>(npart));
+  const auto delta_x = extent_x / npart_T;
+  const auto delta_y = extent_y / npart_T;
+  const auto delta_z = extent_z / npart_T;
   for( auto i = 0ul; i < npart; ++i ) {
-    x_part[i] = bbox_lo[0] + i * delta_x;
-    y_part[i] = bbox_lo[1] + i * delta_y;
-    z_part[i] = bbox_lo[2] + i * delta_z;
+    const auto i_T = traits::from_integer(static_cast<ixx_int>(i));
+    x_part[i] = bbox_lo[0] + i_T * delta_x;
+    y_part[i] = bbox_lo[1] + i_T * delta_y;
+    z_part[i] = bbox_lo[2] + i_T * delta_z;
   }
   x_part.back() = bbox_up[0];
   y_part.back() = bbox_up[1];

include/integratorxx/config.hpp

@@ -0,0 +1,92 @@
+#pragma once
+
+#include <cstddef>
+
+#ifdef ENABLE_STRING_REALS
+#include <string_view>
+#endif
+
+/** @file config.hpp
+ *
+ *  Build-time policy for how numeric literals are spelled and carried.
+ *
+ *  Quadrature code distinguishes two kinds of literal:
+ *
+ *  - **Integral** values, which are exactly representable in every arithmetic
+ *    type the library supports. These are written `IXX_INT(2)` and converted
+ *    with `fp_traits<T>::from_integer`. Note that a value is integral by its
+ *    *value*, not its spelling: a `2.0` in the source is an integral literal.
+ *
+ *  - **Non-integral** values, which are not exactly representable. These are
+ *    written `IXX_REAL(0.6931471805599453) `and converted with
+ *    `fp_traits<T>::from_real`.
+ *
+ *  Keeping the two apart lets arithmetic be carried in exact integers for as
+ *  long as possible, converting once at the end, and lets a ratio of two
+ *  integers go through `fp_traits<T>::divide_integer`, which is a single
+ *  correctly-rounded operation rather than a pre-rounded constant.
+ *
+ *  When `ENABLE_STRING_REALS` is defined, `IXX_REAL` captures the *source text*
+ *  of its argument instead of a `double`. A `double` literal has already lost
+ *  precision by the time any code can inspect it, so a type more precise than
+ *  `double` -- or one that must bound its own error -- cannot recover the
+ *  intended value from it. Carrying the decimal text lets such a type parse it
+ *  directly. See `fp_traits::from_real`.
+ */
+
+namespace IntegratorXX {
+
+/** @brief Integer type for counts, indices, and integral literals.
+ *
+ *  Signed on purpose. Intermediate expressions routinely go negative
+ *  (`2*i - 1`, `npts - 1`, `M - 1`), and an unsigned type would wrap silently
+ *  rather than produce a negative value.
+ */
+using ixx_int = std::ptrdiff_t;
+
+#ifdef ENABLE_STRING_REALS
+/** @brief Non-integral literals carried as their exact decimal source text.
+ *
+ *  `std::string_view` rather than `std::string`: the solid-angle tables are
+ *  `static constexpr`, which requires a literal type.
+ */
+using ixx_real = std::string_view;
+#else
+/// Non-integral literals carried as `double`.
+using ixx_real = double;
+#endif
+
+} // namespace IntegratorXX
+
+/** @brief Denote an integral literal or an integral expression.
+ *
+ *  Parenthesised because expressions are passed, e.g. `IXX_INT(2*i - 1)`.
+ *
+ *  @warning The cast applies to the *result*, not the operands. If the operands
+ *  are unsigned, `IXX_INT(2*i - 1)` evaluates in unsigned arithmetic and only
+ *  then converts, so an underflow has already happened. Convert the count once
+ *  at the top of a scope (`const ixx_int n = IXX_INT(npts);`) and derive the
+ *  rest from it.
+ */
+#define IXX_INT(x) (static_cast<::IntegratorXX::ixx_int>(x))
+
+#ifdef ENABLE_STRING_REALS
+#  define IXX_REAL(x) (::IntegratorXX::ixx_real(#x))
+#else
+/** @brief Denote a non-integral literal.
+ *
+ *  Deliberately expands to its argument unchanged rather than to `(x)`. The
+ *  argument is by contract a single numeric literal token, so no parentheses
+ *  are needed for precedence, and the omission makes the default build's token
+ *  stream byte-identical to one written with bare literals -- which matters
+ *  when the expansion happens ~1.9 million times across the tables.
+ */
+#  define IXX_REAL(x) x
+#endif
+
+namespace IntegratorXX {
+
+/// pi. Spelled as digits because `IXX_REAL(M_PI)` would stringify to "M_PI".
+inline constexpr ixx_real ixx_pi = IXX_REAL(3.14159265358979323846);
+
+} // namespace IntegratorXX

include/integratorxx/generators/impl/impl.hpp

@@ -3,3 +3,11 @@
 #include <integratorxx/generators/impl/radial_factory.hpp>
 #include <integratorxx/generators/impl/s2_factory.hpp>
 
+// Runtime pieces of the molecular_grid API (see
+// integratorxx/molecular_grid/), aggregated here alongside the rest of the
+// non-template runtime generator layer so that a single
+// `#include <integratorxx/generators/impl/impl.hpp>` (exactly once per
+// project) is enough for a full header-only build.
+#include <integratorxx/molecular_grid/impl/defaults.hpp>
+#include <integratorxx/molecular_grid/impl/partition_weights.hpp>
+

include/integratorxx/generators/impl/pruned_grid.hpp

@@ -12,22 +12,23 @@ namespace IntegratorXX {
 
 namespace detail {
 
-template <typename AngularQuadType, typename RadialQuadType>
-auto make_pruned_grid_impl(const RadialQuadType& rq, 
+template <template<typename> class AngularQuadType, typename RadialQuadType>
+auto make_pruned_grid_impl(const RadialQuadType& rq,
   const std::vector<PruningRegion>& pruning_regions) {
 
-  RadialGridPartition<AngularQuadType> rgp;
+  using T = typename RadialQuadType::point_type;
+  RadialGridPartition<AngularQuadType<T>> rgp;
   for( auto& region : pruning_regions ) {
-    rgp.add_quad( rq, region.idx_st, AngularQuadType(region.angular_size) );
+    rgp.add_quad( rq, region.idx_st, AngularQuadType<T>(region.angular_size) );
   }
   rgp.finalize(rq);
 
-  return SphericalGridFactory::generate_pruned_grid(rq, std::move(rgp));
+  return SphericalGridFactory<T>::generate_pruned_grid(rq, std::move(rgp));
 
 }
 
 template <typename RadialQuadType>
-auto make_pruned_grid(const RadialQuadType& rq, 
+auto make_pruned_grid(const RadialQuadType& rq,
   const std::vector<PruningRegion>& pruning_regions) {
 
   if(pruning_regions.size() == 0)
@@ -41,38 +42,39 @@ auto make_pruned_grid(const RadialQuadType& rq,
 
   switch(angular_quad) {
     case AngularQuad::AhrensBeylkin:
-      return make_pruned_grid_impl<ah_type>(rq, pruning_regions);
+      return make_pruned_grid_impl<detail::ah_type>(rq, pruning_regions);
     case AngularQuad::Delley:
-      return make_pruned_grid_impl<de_type>(rq, pruning_regions);
+      return make_pruned_grid_impl<detail::de_type>(rq, pruning_regions);
     case AngularQuad::LebedevLaikov:
-      return make_pruned_grid_impl<ll_type>(rq, pruning_regions);
+      return make_pruned_grid_impl<detail::ll_type>(rq, pruning_regions);
     case AngularQuad::Womersley:
-      return make_pruned_grid_impl<wo_type>(rq, pruning_regions);
+      return make_pruned_grid_impl<detail::wo_type>(rq, pruning_regions);
     default:
       throw std::runtime_error("Unsupported Angular Quadrature");
       abort();
   }
-  
+
 
 }
 
 } // Implementation Details
 
-SphericalGridFactory::spherical_grid_ptr 
-  SphericalGridFactory::generate_pruned_grid( RadialQuad rq, 
-  const RadialTraits& traits, 
-  const std::vector<PruningRegion>& pruning_regions) { 
+template <typename T>
+typename SphericalGridFactory<T>::spherical_grid_ptr
+  SphericalGridFactory<T>::generate_pruned_grid( RadialQuad rq,
+  const RadialTraits<T>& traits,
+  const std::vector<PruningRegion>& pruning_regions) {
 
   switch( rq ) {
 
     case RadialQuad::Becke:
-      return detail::make_pruned_grid( bk_type(traits), pruning_regions );
+      return detail::make_pruned_grid( detail::bk_type<T>(traits), pruning_regions );
     case RadialQuad::MuraKnowles:
-      return detail::make_pruned_grid( mk_type(traits), pruning_regions );
+      return detail::make_pruned_grid( detail::mk_type<T>(traits), pruning_regions );
     case RadialQuad::MurrayHandyLaming:
-      return detail::make_pruned_grid( mhl_type(traits), pruning_regions );
+      return detail::make_pruned_grid( detail::mhl_type<T>(traits), pruning_regions );
     case RadialQuad::TreutlerAhlrichs:
-      return detail::make_pruned_grid( ta_type(traits), pruning_regions );
+      return detail::make_pruned_grid( detail::ta_type<T>(traits), pruning_regions );
 
     default:
       throw std::runtime_error("Unsupported Radial Quadrature");
@@ -83,8 +85,9 @@ SphericalGridFactory::spherical_grid_ptr
 }
 
 
-PrunedSphericalGridSpecification create_pruned_spec(
-  PruningScheme scheme, UnprunedSphericalGridSpecification unp
+template <typename T>
+PrunedSphericalGridSpecification<T> create_pruned_spec(
+  PruningScheme scheme, UnprunedSphericalGridSpecification<T> unp
 ) {
 
   if(!unp.radial_traits) throw std::runtime_error("RadialTraits Not Set");
@@ -93,14 +96,14 @@ PrunedSphericalGridSpecification create_pruned_spec(
       return robust_psi4_pruning_scheme(unp);
     case PruningScheme::Treutler:
       return treutler_pruning_scheme(unp);
-    
+
     // Default to Unpruned Grid
     case PruningScheme::Unpruned:
     default:
       std::vector<PruningRegion> pruning_regions = {
         {0ul, unp.radial_traits->npts(), unp.angular_quad, unp.angular_size}
       };
-      return PrunedSphericalGridSpecification(
+      return PrunedSphericalGridSpecification<T>(
         unp.radial_quad, unp.radial_traits->clone(), pruning_regions
       );
   }