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. */
54export function planeResiduals(x: ArrayLike<number>, y: ArrayLike<number>, z: ArrayLike<number>): Float64Array {
55 const n = x.length;
56 const AtA = new Array(9).fill(0), Atb = [0, 0, 0];
57 for (let i = 0; i < n; i++) {
58 const row = [x[i], y[i], 1];
59 for (let a = 0; a < 3; a++) { Atb[a] += row[a] * z[i]; for (let b = 0; b < 3; b++) AtA[3 * a + b] += row[a] * row[b]; }
60 }
61 const coef = solve3(AtA, Atb);
62 const res = new Float64Array(n);
63 for (let i = 0; i < n; i++) res[i] = z[i] - (coef[0] * x[i] + coef[1] * y[i] + coef[2]);
64 return res;
65}