Revise and document the one-dimensional quadrature code

[garth-wells]

Closes #570.

The helpers that build the one-dimensional rules are documented, and the functions that duplicated each other are removed. The points and weights are unchanged.

Duplication removed

• compute_jacobi_deriv evaluated (P_n^{a,0}) and its derivatives at an array of points, allocating (nderiv + 1, n + 1, npoints) and (n + 1, npoints) arrays. Its only caller asked for one derivative and read a single row of the result. jacobi_value_deriv already computes exactly that from the same recurrence, keeping only the last two values of each sequence, so compute_jacobi_deriv is deleted.
• The map from [-1, 1] to [0, 1] existed in three places: make_quadrature_line, quadrature::gauss_jacobi_rule and quadrature::get_gl_points, each written slightly differently. They now share one helper, gauss_jacobi_rule_01.
• compute_gll_rule is folded into make_gll_line, its only caller.
• The tensor-product loops for quadrilaterals and hexahedra were four copies of the same nested loops, in the Gauss-Jacobi and GLL paths. They now call a single tensor_product helper.

Also dropped: the now-unused mdarray_t alias, an unused local in compute_jacobi_deriv, a commented-out line in compute_gll_rule, and two returns by copy that should have been moves. M_PI, which is POSIX rather than standard C++, is replaced by std::numbers::pi.

Documentation

Every helper now states the domain it works on, the weight function, how many points it uses and the degree it is exact to, and cites the algorithm it implements. The collapsed triangle and tetrahedron rules get the Duffy map written out, along with its Jacobian determinant and the reason the Gauss-Jacobi weight exponents are 1 and 2. The @todo on gauss is resolved. On the public side, get_gll_points now documents the ordering of the points it returns, which was not recorded anywhere: the two endpoints come first, in the order 0, 1, followed by the interior points in increasing order.

The Newton iteration in compute_gauss_jacobi_points keeps its tolerance of 1e-8, now with a comment explaining it: quadratic convergence means an update that small leaves a root accurate to roughly 1e-16. The comment also records that in single precision the update can stagnate at the level of roundoff without reaching the tolerance, so the iteration runs to max_iter with the root already accurate to single precision.

The numbers are unchanged

Checked by dumping every rule from the library built before and after, as hex float:

• all seven cell types × all seven quadrature types × both polyset types × degrees 0 to 24,
• get_gl_points, get_gll_points and gauss_jacobi_rule with a from 0 to 3,
• in both float and double.

That is 3,142,000 double values and 3,133,542 float values, every one of them bit-identical.

Several changes had to be written carefully to keep that true, because the literals in an expression decide whether it is evaluated in double: the original code mixes 1.0 and 0.5 literals into float expressions, so the intermediate arithmetic is done in double and only the result is narrowed. The literals are therefore left as they were.

Tests

Two cases added to cpp/test/quadrature.cpp, covering the two paths that changed shape rather than just gaining comments:

• the quadrilateral and hexahedron rules are the tensor product of the interval rule, with the multi-index ordering the new helper promises;
• get_gl_points(m) agrees with the interval rule that uses m points.

🤖 Generated with Claude Code

[garth-wells]

Pushed a second commit expanding the mathematics in the documentation, and checking it.

Each routine now states the formula it implements rather than only naming the algorithm, so the code can be read against the documentation line by line: the closed forms of the recurrence coefficients and the zeroth moment that fixes beta_0; the Jacobi matrix and weight formula of the Golub-Welsch step; the continued fraction and coefficient replacement in the Lobatto step; the three-term recurrence and its differentiated form, with the local variables renamed to match the symbols; the derivation of the deflated Newton update; the identity and degree of exactness of each rule; the change of variables and degree count for the collapsed triangle and tetrahedron rules.

The statements were checked, not just written. In 50-digit arithmetic: the recurrence against mpmath's Jacobi polynomials; the quadrature weight formula w_i = 2^(a+1) / ((1 - x_i^2) [P_m^(a,0)'(x_i)]^2) for a = 0, 1, 2, 3, 0.5 and 2.75 against the exact integrals; and every stated degree of exactness, confirming in each case that it is attained and that the next degree is not.

That check turned up a real off-by-one in lobatto. The continued fraction for the ratio of consecutive orthogonal polynomials needs alpha[i - 1], but the code used alpha[i] (and alpha[n - 1] instead of alpha[n - 2] in the final step). Against Gautschi's reference implementation:

weight	n	as implemented	fixed	reference
(1-x)^0 (1+x)^0	5	[-1, 1]	[-1, 1]	[-1, 1]
(1-x)^1 (1+x)^0	5	[-1.020924, 0.98995258]	[-1, 1]	[-1, 1]
(1-x)^2 (1+x)^0	5	[-1.0346843, 0.98487798]	[-1, 1]	[-1, 1]
(1-x)^0.5 (1+x)^1.5	5	[-0.99043547, 1.0152268]	[-1, 1]	[-1, 1]

So the prescribed points were not reproduced for any weight that is not symmetric about the origin. The only caller passes the Legendre coefficients, where every alpha is zero and the index therefore does not matter, so this was latent rather than live, and the fix changes no output: the points and weights are still bit-identical to main across the full sweep described above, in both precisions.

With the fix the continued fraction equals -pi_{n-2}(y) / pi_{n-1}(y) for every weight tested, which is what the documentation now claims.

Doxygen builds clean with WARN_AS_ERROR = YES.

[jhale]

This makes the rule useful...

[jhale]

It - Duffy transformed rules therefore require more points than a rule designed specifically for the triangle.

[jhale]

Add: Collapsed rules therefore require more points than a rule designed specifically for the tetrahedron.