Restore the order-39 Ahrens-Beylkin rule: it is 572 points, not 552

[susilehtola]

Fixes #108. Stacked on #101, which withdraws the broken 552-point table;
this restores order 39 with the correct rule.

The 552-point table is missing an orbit

Its weights sum to 0.9632 * 4pi, so it does not integrate a constant
correctly, let alone anything else. Renormalising does not save it — degree 6
is still 0.75% off and degree 8 is 1.4% off, on a rule advertised as order 39.

The exactness conditions are linear in the weights, so the best possible
weights for a fixed point set follow from one least-squares solve. With AB-612
as a control:

grid	best possible weights, full order-39/41 conditions
AB-612 (known good)	8.3e-14, all positive
AB-552	2.4e-02, five negative weights

So no weighting of those 552 points works: the point set itself is incomplete.

The 552 points form ten icosahedral orbits (one of 12, nine of 60). Dividing
the weight deficit, 0.461897446456399, by the candidate orbit sizes leaves one
plausible answer:

12-point orbit -> weight 0.038491454   (above the existing weight range)
20-point orbit -> weight 0.023094872   <-- sits at the top of it
30-point orbit -> weight 0.015396582   (below the range)
60-point orbit -> weight 0.007698291   (far below)

Adding the 20-point dodecahedral-vertex orbit gives 572 = 12 + 20 + 9*60,
and the same least-squares then reproduces the file's own ten weights exactly
while assigning the new orbit precisely deficit / 20: 3.674e-13, all weights
positive. That is the missing orbit, not a coincidental fit.

Verification

Refined against the full order-39 condition set:

1.529e-15 -> 8.862e-28 -> 7.905e-52     worst error 3.338e-51

Emitted at double precision, in the same format as the other tables, it
integrates as well as its neighbours:

AB-492 shipped (control)   order 37   8.331e-16
AB-572 new                 order 39   7.605e-16
AB-612 shipped (control)   order 41   9.003e-16

npts_by_algebraic_order(39) returns 572, algebraic_order_by_npts(572)
returns 39, next_algebraic_order no longer skips 39, and asking for 552 still
throws. The polynomial-exactness test added in #101 covers the new size
automatically. Library suite 9/9.

Provenance

The original file received from Beylkin, qsph1-39-552DP.dat, has the same
defect: sum(w)/4pi = 0.963243369097, 5.7e-02 relative error at degree 8. So
IntegratorXX's conversion was faithful and the orbit was already missing
upstream — worth reporting back. Every other one of the 79 files in that set is
exact to ~1e-16; order 39 is the only bad one.

Confirmed against the paper

Ahrens & Beylkin, Rotationally invariant quadratures for the sphere,
Proc. R. Soc. A 465, 3103 (2009), figure 3 caption:

there are nine generators with orbits of length 60 and one node at the face
centre with orbit length 20, yielding 28 nodes in the interior of the
triangle and 572 total nodes on the sphere

which is exactly the reconstructed structure: 9x60 + 20 (face centre) + 12
(icosahedron vertices) = 572. So this configuration is a documented member of
the family, not an exception — my earlier concern that 572 = 12 + 20 + 9*60
broke the 12 + k*60 pattern of the other 56 tabulated sizes was unfounded;
that pattern just reflects which templates happen to be tabulated.

The algebraic order is 39 exactly, matching the filename. Measured on the
refined rule at 60 digits:

degree   worst |quad - exact|
   36        1.278e-50   exact
   37        2.102e-52   exact
   38        1.310e-50   exact
   39        2.116e-52   exact
   40        7.028e-13   fails
   41        4.155e-14   fails

🤖 Generated with Claude Code

[wavefunction91]

Two things on the commit message, since that's what lands in history:

• It ends by leaving open whether this is the published rule. The description settles that, and so does the paper's Eq. (3.2): the Qvf quadrature has the 12 vertices, one fixed generator at the face centres (20 nodes), and 60-point orbits, "with a total of 60Ng + 32 nodes", which is 572 for Ng = 9. Could the message say so rather than end on the open question?
• The Fig. 3.3 caption quoted in the description continues: "The resulting quadrature found using Newton's method integrates the subspace of maximum order and degree N = 40." This table is exact through 39 (relative residual ≤ 2.1e-16) and fails at 40 (8e-3). That fits Beylkin's file name and the paper's p. 9 note that, when the full system appears to have no solution, they "reduce the number of equations by removing those from the subspace of the highest degree". So 39 looks like the right label, but a reader who follows the citation will find 40, and a sentence on it would save them the detour.

(Quotes are from the authors' posted copy, https://amath.colorado.edu/faculty/beylkin/papers/AHR-BEY-2009net.pdf, where the figure is numbered 3.3.)

One more on the message: the "2.4e-02, five negative weights" figure for the 552 point set depends on how the conditions are scaled, which isn't stated. With orthonormal harmonics through degree 39 I get a residual of 1.10 and no negative weights; with unnormalised monomials 3.1e-05 and two negative. AB-612 sits near 1e-14 in every basis, so the conclusion (no exact weighting exists) holds either way; it's just the specific numbers that need their basis.

[wavefunction91]

This is the only table in the directory that isn't a transcription of one of Beylkin's qsph1 files, so a line here would help a future reader who opens the file rather than the history: reconstructed from qsph1-39-552DP.dat by adding the face-centre orbit, then refined. It's the paper's Qvf construction (Ahrens & Beylkin 2009, Eq. 3.2: "60Ng + 32 nodes", with Ng = 9).