Repository navigation
Ahrens-Beylkin 552-point grid is not a valid order-39 rule #108
Description
The tabulated 552-point Ahrens-Beylkin grid
(include/integratorxx/quadratures/s2/ahrens_beylkin/ahrens_beylkin_552.hpp)
is not a valid quadrature rule. It is advertised as algebraic order 39, but:
- its weights sum to
12.104473167902746, i.e.0.9632 * 4pi, so it does
not even integrate a constant correctly; - after renormalising to
4piit integrates monomials up to degree 4
exactly, but is still wrong at degree 6 (7.5e-03relative) and degree 8
(1.4e-02relative).
So it is not merely mis-normalised -- renormalising would produce a grid
that looks correct at low order and is silently wrong at high order, which
is worse than the current state.
Every other size in the family is exact to machine precision:
AB 432 sum(w)=12.566370614359144 sum/4pi=1.000000000
AB 492 sum(w)=12.566370614359053 sum/4pi=1.000000000
AB 552 sum(w)=12.104473167902746 sum/4pi=0.963243369 <---
AB 612 sum(w)=12.566370614359178 sum/4pi=1.000000000
AB 672 sum(w)=12.566370614359140 sum/4pi=1.000000000
I swept all 230 grid sizes advertised by the npts_by_algebraic_order
tables across all four angular families, checking exact integration of even
monomials up to each grid's own claimed algebraic order (capped at degree 8
for run time). This is the only failure:
Lebedev : 32 sizes checked, 0 bad
Delley : 18 sizes checked, 0 bad
AhrensBeylkin : 55 sizes checked, 1 bad <-- 552
Womersley : 125 sizes checked, 0 bad
Reproduction:
#include <integratorxx/quadratures/s2/ahrens_beylkin.hpp>
#include <numeric>
#include <cstdio>
int main() {
IntegratorXX::AhrensBeylkin<double> q(552);
auto& w = q.weights();
printf("sum(w) = %.15f (4pi = %.15f)\n",
std::accumulate(w.begin(), w.end(), 0.0), 4*M_PI);
}All 552 points are distinct and of unit norm, and the weights are
individually of plausible magnitude, so the table looks like it lost or
duplicated a subset of the rule rather than being corrupted wholesale. I
have not been able to determine the exact provenance of the error.
Interim action taken. PR #101 withdraws the size from the
dispatch and from the algebraic order tables, so requesting 552 points now
throws rather than returning a 3.7%-wrong quadrature, and
next_algebraic_order rounds order 39 up to 41. The data file is left in
place.
What is needed. Regenerating the table from the original Ahrens-Beylkin
data (Proc. R. Soc. A 465, 3103 (2009)) and re-enabling the size. The
polynomial-exactness test added in the same PR will confirm a correct
replacement.
Activity
susilehtola commentedon Sep 1, 2026
Resolved by #113: the order-39 rule is 572 points, not 552, and the table
was missing its 20-point (dodecahedral vertex) orbit.
The exactness conditions are linear in the weights, so the best possible
weights for a fixed point set follow from a single least-squares solve. With
AB-612 as a control:
No weighting of the 552 points works. Dividing the weight deficit by candidate
orbit sizes leaves exactly one plausible answer — a 20-point orbit at weight
0.023094872, right at the top of the existing range. Adding the dodecahedral
orbit gives 572 = 12 + 20 + 9*60, and the same solve then reproduces the file's
own ten weights exactly while assigning the new orbit precisely deficit/20.
Refined against the full order-39 set: 1.529e-15 -> 8.862e-28 -> 7.905e-52.
Provenance correction. I suggested above that this looked like a
transcription fault. It is not: the original file received from Beylkin,
qsph1-39-552DP.dat, has the same defect (sum(w)/4pi = 0.963243369097,5.7e-02 at degree 8). IntegratorXX's conversion was faithful and the orbit was
already missing upstream. Every other one of the 79 files in that set is exact
to ~1e-16, so this is worth reporting back to the authors.
One open question is carried into #113: all 56 other sizes in the family are of
the form
12 + k*60, so 572 is structurally an exception and the naming shouldbe checked against the paper. The rule itself is verified regardless.
susilehtola commentedon Sep 1, 2026
Confirmed against the paper, so the open question in my previous comment is
closed.
Ahrens & Beylkin, Proc. R. Soc. A 465, 3103 (2009), figure 3 caption:
That is exactly the reconstructed structure — 9x60 + 20 (face centre) + 12
(icosahedron vertices) = 572. My earlier worry that 572 broke the family's
12 + k*60pattern was unfounded: face-centred templates are a documentedconstruction, and the pattern across the other 56 tabulated sizes just reflects
which templates were tabulated.
The order is 39 exactly, matching the filename. On the refined rule at 60
digits it is exact through degree 39 (2.1e-52) and fails at 40 (7.0e-13).
So
qsph1-39-552DP.datis an order-39, 572-node, face-centred rule that lostits 20-node face-centre orbit; the filename records the 552 rows that survived.