Function gauss_legendre 

pub fn gauss_legendre<T: Real>(n: usize) -> Result<GaussLegendre<T>>

Computes the n-point Gauss–Legendre rule, exact for polynomials of degree ≤ 2n − 1.

Nodes ascend, and the rule is exactly symmetric: nodes[i] == -nodes[n - 1 - i] and weights[i] == weights[n - 1 - i]; for odd n the middle node is +0. Two methods:

• Newton (n ≤ 100, and every n for DoubleDouble precision). The rule is computed in double-double and converted with Real::from_f64_pair, so it is rounded once: the f64 rule is the double-double rule rounded to nearest.
• Asymptotic (n > 100, for f64 precision: f64 and Dual<f64, N>). Bogaert’s (2014) expansions of the angle θₖ of each node and of its weight, in powers of 1/(n + ½)², around the zeros of the Bessel function J₀; the weights come from the angle, not from the node, so they keep the accuracy of the nodes next to ±1.

Dual rules are constants (zero derivatives) whose values are bitwise those of f64 or DoubleDouble. No platform libm is used (the starting guesses and the expansions use the double-double sine and cosine), only IEEE 754 operations and fused multiply-adds, so the bits of the rule do not depend on the platform.

Accuracy

Verified against Arb for every node of the rules with n ≤ 1000, and for the nodes nearest to ±1 and to 0 at n = 5000, 20 000 and 2¹⁷ (the cap), and against the double-double rule for every node of every rule with 101 ≤ n ≤ 3000 (an ignored test, which CI runs weekly) and of eleven rules up to n = 4001:

• f64, n ≤ 100: node absolute error ≤ 2⁻⁵², weight relative error ≤ 2⁻⁵². Measured: at most 0.25·2⁻⁵² and 0.46·2⁻⁵², the rounding of the double-double rule.
• f64, n > 100: node absolute error ≤ 32·2⁻⁵², weight relative error ≤ 32·2⁻⁵². Measured: at most 2.8·2⁻⁵² and 4.0·2⁻⁵².
• DoubleDouble: node absolute error ≤ 2⁻⁹⁸, weight relative error ≤ 2⁻¹⁰⁰·(4 + n²). The n² comes from the conditioning of 1 − x² next to ±1: a node error δ changes the weight by 2|x|δ/(1 − x²) relative, and 1 − x² ≈ 5.8/n² at the extreme nodes. Measured: at most 0.51·2⁻¹⁰⁶ for nodes; for weights about 2⁻⁸⁵, 2⁻⁸² and 2⁻⁸⁰ at n = 5000, 20 000 and 2¹⁷.

Cost

O(n) memory. Asymptotic: O(n) time, about 0.3 µs per point, measured on an Apple M4 Pro as 0.3 ms at n = 1000, 6 ms at n = 20 000 and 40 ms at the cap. Newton: O(n²) time; see MAX_GAUSS_LEGENDRE_POINTS for the cost model. Each of the ⌈n/2⌉ roots takes four to seven evaluations of the three-term recurrence in double-double (about four on average for n ≥ 5000). Measured on an Apple M4 Pro: 0.8 ms at n = 100, 32 ms at n = 1000, 0.72 s at n = 5000 and 11 s at n = 20 000. The call cannot be interrupted.

Errors

• Error::InvalidInput if n == 0 or n > MAX_GAUSS_LEGENDRE_POINTS.
• Error::NotConverged if a Newton iteration fails to converge, with achieved the last Newton step |Δx| (absolute, in x; +∞ if it was NaN) and requested 2⁻⁹⁰. The tests have not found any n where this happens: they check every root of every rule with n ≤ 300 (and, in an ignored test, n ≤ 3000), and the roots nearest ±1 and 0 for over a hundred values of n up to the cap.

Examples

use scatrs::special::gauss_legendre;

// The 3-point rule integrates x⁴ exactly: ∫ x⁴ dx over [-1, 1] = 2/5.
let rule = gauss_legendre::<f64>(3)?;
let integral: f64 = rule.nodes.iter().zip(&rule.weights).map(|(x, w)| w * x.powi(4)).sum();
assert!((integral - 0.4).abs() < 1e-15);