1/** Gauss-Legendre nodes (ascending) and weights on [-1, 1], by Newton iteration. */
2export function leggauss(n: number): { x: Float64Array; w: Float64Array } {
3 const x = new Float64Array(n);
4 const w = new Float64Array(n);
5 const m = (n + 1) >> 1;
6 for (let i = 0; i < m; i++) {
7 let z = Math.cos((Math.PI * (i + 0.75)) / (n + 0.5));
8 let pp = 0;
9 for (let iter = 0; iter < 100; iter++) {
10 let p1 = 1;
11 let p2 = 0;
12 for (let j = 1; j <= n; j++) {
13 const p3 = p2;
14 p2 = p1;
15 p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
16 }
17 pp = (n * (z * p1 - p2)) / (z * z - 1);
18 const z1 = z;
19 z = z1 - p1 / pp;
20 if (Math.abs(z - z1) < 1e-15) break;
21 }
22 x[i] = -z;
23 x[n - 1 - i] = z;
24 w[i] = w[n - 1 - i] = 2 / ((1 - z * z) * pp * pp);
25 }
26 return { x, w };
27}
29/** P_0(x), ..., P_norder(x) for each x: array of shape (norder+1, n), index k*n + i. */
30export function legendreVander(x: ArrayLike<number>, norder: number): Float64Array {
31 const n = x.length;
32 const P = new Float64Array((norder + 1) * n);
33 for (let i = 0; i < n; i++) {
34 P[i] = 1;
35 if (norder >= 1) P[n + i] = x[i];
36 for (let k = 1; k < norder; k++) {
37 P[(k + 1) * n + i] = ((2 * k + 1) * x[i] * P[k * n + i] - k * P[(k - 1) * n + i]) / (k + 1);
38 }
39 }
40 return P;
41}