Re-express the Chebyshev weights to avoid cancellation

[susilehtola]

Refs #118. Stacked on #119, which it needs only to avoid a textual conflict in
the same files.

Three of the five sites from the #118 audit — the ones with exact analytic
re-expressions and no design decision attached.

All three Chebyshev rules form a quantity that cancels as the nodes approach
the interval ends, and in each case the node is a cosine, so the cancelling
quantity is available exactly:

rule	as written	exactly
GaussChebyshev1	(pi/n) sqrt(1-x^2), x = cos t	(pi/n) sin t
GaussChebyshev2	h sqrt(1-x^2), x = cos t	h sin t
GaussChebyshev3	2h x sqrt((1-x)/x), x = cos^2 t	h sin(2t)

The third collapses further than the first two: 2h cos^2(t) tan(t) is
2h cos(t) sin(t), so the whole product is h sin(2t).

Measured

Worst-case relative error in the weights, against a 50-digit reference:

rule	n=50	n=200	n=800
GaussChebyshev1	5.55e-14 → 2.15e-15	9.44e-13 → 4.65e-14	1.78e-11 → 1.13e-13
GaussChebyshev2	5.87e-15 → 1.67e-15	9.02e-14 → 1.05e-14	2.93e-12 → 1.34e-13
GaussChebyshev3	2.44e-13 → 6.86e-15	6.74e-13 → 2.08e-14	2.81e-11 → 4.33e-14

The error grows as O(n^2) before the change, because 1 - cos^2 t is O(1/n^2)
against terms of O(1) — so it is worst exactly where a large radial grid is most
likely to be trusted and least likely to be checked.

Two caveats about what this does and does not buy

It changes no integral at double precision. These rules have their weights
rescaled by sqrt(1-x^2), so they are not Gauss rules for polynomials and their
truncation error dominates. Integrating 1 and x^2 gives identical results
before and after, at 1e-4 to 1e-7. I checked against master specifically to see
whether the fix was observable at that level, and it is not. The gain is the
fidelity of the weights to their defining formula, which matters to a consumer
that amplifies them — a radial Jacobian of 1/(1-x) or 1/(1-x)^2 — and to
anything comparing against a reference.

It is not a double-precision artefact. The loss is relative and therefore
independent of working precision: at 40 digits, GaussChebyshev1 at n=800 still
loses two orders, 5.62e-37 against 4.11e-39 re-expressed. Any
arbitrary-precision use of these rules inherits it, which is the strongest
argument for the change.

Why the tests did not catch it

They check weights with an absolute tolerance of 10*eps, and the weights
near the endpoints are O(1e-4), so the effective relative tolerance is ~2e-11 —
looser than the error being measured. The suite passes either way (9/9), and
the reference values in test/gausschebyshev{1,2}.cxx are unchanged.

Left open deliberately

The other two audit sites each need a decision rather than a local edit:

• The radial transforms' own 1-x (becke, treutlerahlrichs, mhl,
  muraknowles). The transform cannot re-express it — only the primitive knows
  the node is a cosine, where 1 - cos t = 2 sin^2(t/2) is exact. So either the
  primitive supplies the complements with the node, or transforms are written
  against t. Worth 2.68e-12 → 4.58e-17 at n=800, and 5.4e-12 in the Becke
  Jacobian. This is an interface change touching all four transforms, and it
  overlaps the restructuring proposed in Add an arbitrary-precision generator for the angular grid tables #110.
• The GaussChebyshev2Modified node (GaussChebyshev2Modified: the outermost node loses ~4 log10(n) digits, breaking logarithmic radial maps #117). No exact re-expression; it
  needs a series in u and a branch, and it changes node values, so it
  breaks bit-for-bit ERKALE parity for downstream consumers.

🤖 Generated with Claude Code

https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF

[wavefunction91]

See inline comments.

[wavefunction91]

This removes the 1 - x*x cancellation, but the O(n) growth in the "after" column comes from where sin is now evaluated. The loop computes the half with i = npts - idx, i.e. t ∈ (π/2, π), and reflects it. Near π, the absolute rounding of ti (≈ ulp(π)) becomes an O(nε) relative error in sin t ≈ π − t, and the reflection copies it to both endpoints. GC3 has the same issue as 2t → π.

The complementary angle is available exactly from the index: sin(tᵢ) = sin((2·idx+1)·pi_ov_2npts) here, sin((idx+1)·pi_ov_npts_p_1) in GC2 (line 84), and sin(min(2i−1, 2n+2−2i)·pi_ov_2n_p_1) in GC3 (line 85). Against a 50-digit reference, that gives a worst relative weight error of 3.6e-16 / 4.7e-16 / 2.8e-16 at n=800 (vs 1.1e-13 / 1.3e-13 / 4.3e-14 on this head), flat in n, and leaves the nodes untouched.

With that, a relative-tolerance check on the endpoint weights at large n (e.g. 1e-14 at n=800) would pin the improvement; the existing absolute 10*eps checks can't see it, as the description notes.

[susilehtola]

You are right, and thank you for the diagnosis — I had attributed all of the remaining growth to the 1 - x*x cancellation, when the dominant term was this second, independent one.

Fixed in 6b8c630 using the forms you give. Worst relative weight error against a 50-digit reference at n=800, which reproduces your figures:

rule	before	after
T1	1.1e-13	3.6e-16
T2	1.3e-13	4.7e-16
T3	4.3e-14	2.8e-16

Flat in n now rather than growing (T1 is 1.8e-16 at n=50 and 2.5e-16 at n=100), and the nodes are untouched. For T3 the reflection is written as k = (m <= npts) ? m : 2*npts + 1 - m with m = 2i - 1, since m h passes pi/2 exactly when m > npts.

The new "Gauss-Chebyshev extremal weights at large N" case pins all three at a relative 1e-14 at n=800. Its reference is the Taylor series of sin about zero rather than std::sin, so it does not assume the implementation it is checking; at these angles the first omitted term is O(h^7), around 1e-22 relative. Five of its six assertions fail on the previous head — the one that passes is the T3 node nearest zero, which was already evaluated at a small angle.

Note this is now rebased onto the updated #119 branch, which gained a commit of its own.

[wavefunction91]

See inline comment.

[wavefunction91]

The labels are swapped: front() is the node nearest 0 (x ≈ 3.9e-6) and carries h·sin 2h; back() (x ≈ 1) carries h·sin h. The assertions are right.

[susilehtola]

Correct, and thank you — t = pi/2 - h exactly at idx = 0, so front() is the node nearest 0 with x = sin^2(h), about 3.85e-6 at n=800, matching your figure. Comment fixed in e45f43f; the assertions were already paired the right way.

[wavefunction91]

#119 is in as f4af7ca. Since it was squash-merged, this branch still carries #119's original commits, including the test comment from before ef50142, and that block now conflicts with master. git rebase --onto origin/master 00061bc1 replays just the three commits from c638bca on, and applies cleanly here. Once that's pushed and green, I'll merge.

[susilehtola]

Rebased onto e56a92c now that #119 has merged. #119 was squash-merged, so this branch still carried the two pre-squash commits and had gone conflicting; both are dropped, leaving only the three commits that belong to this PR. The diff against master is unchanged in substance, and master's corrected comments about the MK/MHL aliases are picked up from there rather than duplicated here.

10/10 tests pass on the rebased branch, including the assoc_legendre case that arrived with #112.