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