Repository navigation
Interest in a Floating-Point Agnostic Implementation? #100
Description
Long story short, NWChemEx needs DFT grids that work with arbitrary floating point (FP) types. We'd like to use IntegratorXX, but AFAIK the current blockers are:
-
constexprtables for angular components only work for literal FPs. It's fine to store the tables asdoubleor whatever, but there needs to be some logic that will convert them to the user's specified FP type. Solution would be to run them through some template meta-programming to figure out if we use the table directly or if we need to convert them. -
The use of
std::expand the like. Custom FPs often provide their own versions of the transcendental functions, we'd either need to remove thestd::and rely on argument dependent lookup (assuming the custom FP use the same signature) or write traits classes that can be customized by the user to control what gets called (default implementation uses ADL). IMHO the traits class is the better option. -
There are
std::is_floating_point_vgates on thequadrature_traits<MuraKnowles<PointType, WeightType>that would need to be removed. -
Literals, i.e., expressions like
x + 1.0. I think it's reasonable to assume custom FP will defineoperator+(double), but there may be places where wrapping1.0asT(1.0)is needed (if going with traits for 2, can just make a method for convertingdoubleobjects). -
Somewhat ironically, the Gauss-Legendre Newton iterations need to be pinned to a literal floating-point type like
doubleand the result converted to the custom FP type. The iterations are ultimately an auxiliary calculation to find a root. It's the value of the root, not how it was obtained, that the grid cares about. Thus, for the custom FP types I'm worried about, we don't want the FP type to enter into the iterations.
I need to make these changes for IntegratorXXX to work for us. My question is once I've made these changes: are you interested in having them contributed back?
Activity
susilehtola commentedon Aug 26, 2026
I am personally heading in this direction already in my projects, where the floating point datatype is templated so that I can use e.g. float128 or MPFR datatypes, and I would really like to see this in IntegratorXX as well.
ryanmrichard commentedon Aug 27, 2026
FWIW the changes weren't as bad as I first thought. According to Claude ryanmrichard#1 would satisfy this issue. It needs some more verification/clean-up before it is mergeable, but that gives the IntegratorXX maintainers (@wavefunction91, but not sure if there's others) an idea of what the push to this repo would look like.
wavefunction91 commentedon Oct 13, 2026
Thanks both for pushing on this. Before #110 and #115 go further, here's the direction I'd like for the representation, since it constrains both. IntegratorXX sits under GauXC, and through it qdk-chemistry, so higher precision has to leave the float/double grids bit-for-bit unchanged and can't grow the headers. That rules out storing the tables as strings, and build options that change the type of the shipped tables.
floatanddoubleuse term 0 only, so their grids don't move; wider types sum the terms in T, chosen by an explicit trait rather than by digit count. Two terms can't reach binary128 even in the best case (about 106 bits); three cover it and the 40 digits Add an arbitrary-precision generator for the angular grid tables #110 emits, with the emitter verifying every split. That's a fixed ceiling, not arbitrary precision, and the tool can emit more terms locally for anyone who needs them. At three terms, LL and Delley come to 22,848 doubles.long double/ Boost.Multiprecision test, then Delley the same way. Ahrens-Beylkin stays as it is: its orbits are under the rotation group, and regenerating its points from rotations changes their double values and order. Womersley has no symmetry to exploit and stays double-only. For both, the tool can still produce higher-precision tables locally.I'd rather not refine at runtime: it would put a rank-revealing least-squares solve in T inside the library, against the 2e24 and 8e26 condition numbers #110 reports for delley_1730 and delley_2030.
The per-grid structs (
LebedevLaikovGrids::lebedev_laikov_110<T>and friends) are installed headers, so (1) is a source-API change even thoughLebedevLaikov<T>itself doesn't change. That's fine, but the PR making it should say so.