1// Small dense linear algebra: symmetric eigenvalues by cyclic Jacobi, and a least-squares plane fit.
3/** Eigenvalues (ascending) of a symmetric n x n matrix given row-major; the input is copied. */
4export function symmetricEigenvalues(A: Float64Array | number[], n: number): Float64Array {
5 const a = Float64Array.from(A);
6 for (let sweep = 0; sweep < 60; sweep++) {
7 let off = 0;
8 for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += a[p * n + q] * a[p * n + q];
9 let diag = 0;
10 for (let p = 0; p < n; p++) diag += a[p * n + p] * a[p * n + p];
11 if (off <= 1e-30 * (diag + 1e-300)) break;
12 for (let p = 0; p < n; p++) {
13 for (let q = p + 1; q < n; q++) {
14 const apq = a[p * n + q];
15 if (apq === 0) continue;
16 const app = a[p * n + p], aqq = a[q * n + q];
17 const theta = (aqq - app) / (2 * apq);
18 const t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));
19 const c = 1 / Math.sqrt(t * t + 1), s = t * c;
20 for (let k = 0; k < n; k++) {
21 const akp = a[k * n + p], akq = a[k * n + q];
22 a[k * n + p] = c * akp - s * akq; a[k * n + q] = s * akp + c * akq;
23 }
24 for (let k = 0; k < n; k++) {
25 const apk = a[p * n + k], aqk = a[q * n + k];
26 a[p * n + k] = c * apk - s * aqk; a[q * n + k] = s * apk + c * aqk;
27 }
28 }
29 }
30 }
31 const ev = new Float64Array(n);
32 for (let i = 0; i < n; i++) ev[i] = a[i * n + i];
33 return ev.sort();
34}
36/** Solve the 3 x 3 system A x = b (row-major A) by Gaussian elimination with pivoting. */
37export function solve3(A: number[], b: number[]): number[] {
38 const M = [A.slice(0, 3).concat(b[0]), A.slice(3, 6).concat(b[1]), A.slice(6, 9).concat(b[2])];
39 for (let c = 0; c < 3; c++) {
40 let piv = c;
41 for (let r = c + 1; r < 3; r++) if (Math.abs(M[r][c]) > Math.abs(M[piv][c])) piv = r;
42 [M[c], M[piv]] = [M[piv], M[c]];
43 if (Math.abs(M[c][c]) < 1e-300) throw new Error("singular 3x3 system");
44 for (let r = 0; r < 3; r++) {
45 if (r === c) continue;
46 const f = M[r][c] / M[c][c];
47 for (let k = c; k < 4; k++) M[r][k] -= f * M[c][k];
48 }
49 }
50 return [M[0][3] / M[0][0], M[1][3] / M[1][1], M[2][3] / M[2][2]];
51}
53/** Least-squares fit z = a x + b y + c through points (x_i, y_i, z_i); returns the residuals z - fit.
54 *
55 * The fit is done by Gram-Schmidt on the columns (1, x, y), dropping a column that is (nearly)
56 * dependent on the previous ones, which is what numpy's minimum-norm least squares does in the
57 * Python package: the collinear points of a stick's line get a line fit, a single point a zero
58 * residual. The normal equations would be singular in those cases. */
59export function planeResiduals(x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>): Float64Array {
60 const n = x.length;
61 const res = Float64Array.from(z as number[]);
62 const basis: Float64Array[] = [];
63 for (const col of [Array.from({ length: n }, () => 1), Array.from(x), Array.from(y)]) {
64 const v = Float64Array.from(col);
65 let norm0 = 0;
66 for (let i = 0; i < n; i++) norm0 += v[i] * v[i];
67 for (const b of basis) {
68 let d = 0;
69 for (let i = 0; i < n; i++) d += b[i] * v[i];
70 for (let i = 0; i < n; i++) v[i] -= d * b[i];
71 }
72 let norm = 0;
73 for (let i = 0; i < n; i++) norm += v[i] * v[i];
74 if (norm === 0 || norm <= 1e-20 * norm0) continue;
75 const inv = 1 / Math.sqrt(norm);
76 for (let i = 0; i < n; i++) v[i] *= inv;
77 basis.push(v);
78 }
79 for (const b of basis) {
80 let d = 0;
81 for (let i = 0; i < n; i++) d += b[i] * res[i];
82 for (let i = 0; i < n; i++) res[i] -= d * b[i];
83 }
84 return res;
85}