1/**
2 * Gauss-Legendre quadrature nodes and weights, computed in double
3 * precision by Newton iteration on P_n (cf. gauss_nodes() in SHTNS
4 * sht_legendre.c).
5 *
6 * Returns nodes x_i = cos(theta_i) in DECREASING order (theta increasing,
7 * north pole first), and weights w_i for integration over x in [-1, 1]:
8 * integral f(x) dx ~= sum_i w_i f(x_i), exact for polynomials of
9 * degree <= 2n - 1.
10 */
11export function gaussNodesWeights(n: number): { x: Float64Array; w: Float64Array } {
12 const x = new Float64Array(n);
13 const w = new Float64Array(n);
14 const m = (n + 1) >> 1;
15 for (let i = 0; i < m; i++) {
16 // initial guess (Tricomi-like), then Newton
17 let z = Math.cos((Math.PI * (i + 0.75)) / (n + 0.5));
18 let pp = 0;
19 for (let iter = 0; iter < 100; iter++) {
20 // evaluate P_n(z) and P_{n-1}(z) by recurrence
21 let p1 = 1.0;
22 let p2 = 0.0;
23 for (let j = 1; j <= n; j++) {
24 const p3 = p2;
25 p2 = p1;
26 p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
27 }
28 pp = (n * (z * p1 - p2)) / (z * z - 1.0);
29 const dz = p1 / pp;
30 z -= dz;
31 if (Math.abs(dz) < 1e-15 * Math.abs(z) + 1e-300) {
32 // one extra iteration for full convergence
33 let q1 = 1.0, q2 = 0.0;
34 for (let j = 1; j <= n; j++) {
35 const q3 = q2; q2 = q1;
36 q1 = ((2 * j - 1) * z * q2 - (j - 1) * q3) / j;
37 }
38 pp = (n * (z * q1 - q2)) / (z * z - 1.0);
39 z -= q1 / pp;
40 break;
41 }
42 }
43 x[i] = z; // largest roots first => theta increasing
44 x[n - 1 - i] = -z;
45 const wi = 2.0 / ((1.0 - z * z) * pp * pp);
46 w[i] = wi;
47 w[n - 1 - i] = wi;
48 }
49 if (n & 1) x[m - 1] = 0.0; // exact for odd n
50 return { x, w };
51}