// Dense symmetric eigensolver: Householder reduction to tridiagonal form followed by the
// implicit QL algorithm with accumulation of the transformations (the tred2 / tql2 pair of
// EISPACK, in the form given by Numerical Recipes), rearranged so that every inner loop runs
// along a contiguous row. The eigenvectors are returned as the rows of vectorsT.

const EPS = Number.EPSILON;

function pythag(a: number, b: number): number {
  return Math.hypot(a, b);
}

/** Householder tridiagonalisation of the symmetric row-major n x n matrix a (destroyed). On
 * return d holds the diagonal, e the sub-diagonal (e[0] = 0) and a the orthogonal matrix Q. */
function tred2(a: Float64Array, n: number, d: Float64Array, e: Float64Array, onProgress?: (frac: number) => void): void {
  const g = new Float64Array(n);
  for (let i = n - 1; i >= 1; i--) {
    const l = i - 1;
    const ri = i * n;
    let h = 0, scale = 0;
    if (l > 0) {
      for (let k = 0; k <= l; k++) scale += Math.abs(a[ri + k]);
      if (scale === 0) {
        e[i] = a[ri + l];
      } else {
        for (let k = 0; k <= l; k++) {
          a[ri + k] /= scale;
          h += a[ri + k] * a[ri + k];
        }
        let f = a[ri + l];
        const gg = f >= 0 ? -Math.sqrt(h) : Math.sqrt(h);
        e[i] = scale * gg;
        h -= f * gg;
        a[ri + l] = f - gg;
        f = 0;
        // p = A u / h, using full rows of the (symmetric) leading block
        for (let j = 0; j <= l; j++) {
          a[j * n + i] = a[ri + j] / h;
          let s = 0;
          const rj = j * n;
          for (let k = 0; k <= l; k++) s += a[rj + k] * a[ri + k];
          e[j] = s / h;
          f += e[j] * a[ri + j];
        }
        const hh = f / (h + h);
        for (let j = 0; j <= l; j++) e[j] -= hh * a[ri + j];
        // rank-two update of the full leading block: A -= u q^T + q u^T
        for (let j = 0; j <= l; j++) {
          const fj = a[ri + j], gj = e[j];
          const rj = j * n;
          for (let k = 0; k <= l; k++) a[rj + k] -= fj * e[k] + gj * a[ri + k];
        }
      }
    } else {
      e[i] = a[ri + l];
    }
    d[i] = h;
    if (onProgress && (i & 63) === 0) onProgress((n - i) / n);
  }
  d[0] = 0;
  e[0] = 0;
  // accumulate the transformations into a (Q)
  for (let i = 0; i < n; i++) {
    const l = i - 1;
    const ri = i * n;
    if (d[i] !== 0) {
      g.fill(0, 0, l + 1);
      for (let k = 0; k <= l; k++) {
        const aik = a[ri + k];
        if (aik === 0) continue;
        const rk = k * n;
        for (let j = 0; j <= l; j++) g[j] += aik * a[rk + j];
      }
      for (let k = 0; k <= l; k++) {
        const aki = a[k * n + i];
        if (aki === 0) continue;
        const rk = k * n;
        for (let j = 0; j <= l; j++) a[rk + j] -= g[j] * aki;
      }
    }
    d[i] = a[ri + i];
    a[ri + i] = 1;
    for (let j = 0; j <= l; j++) {
      a[j * n + i] = 0;
      a[ri + j] = 0;
    }
  }
}

/** Implicit QL on the tridiagonal (d, e), with the rotations applied to the rows of zT. */
function tql2(d: Float64Array, e: Float64Array, n: number, zT: Float64Array, onProgress?: (frac: number) => void): void {
  for (let i = 1; i < n; i++) e[i - 1] = e[i];
  e[n - 1] = 0;
  for (let l = 0; l < n; l++) {
    let iter = 0;
    let m: number;
    do {
      for (m = l; m < n - 1; m++) {
        const dd = Math.abs(d[m]) + Math.abs(d[m + 1]);
        if (Math.abs(e[m]) <= EPS * dd) break;
      }
      if (m !== l) {
        if (iter++ === 100) throw new Error("tql2: too many iterations");
        let g = (d[l + 1] - d[l]) / (2 * e[l]);
        let r = pythag(g, 1);
        g = d[m] - d[l] + e[l] / (g + (g >= 0 ? Math.abs(r) : -Math.abs(r)));
        let s = 1, c = 1, p = 0;
        let i: number;
        for (i = m - 1; i >= l; i--) {
          const f = s * e[i], b = c * e[i];
          e[i + 1] = r = pythag(f, g);
          if (r === 0) {
            d[i + 1] -= p;
            e[m] = 0;
            break;
          }
          s = f / r;
          c = g / r;
          g = d[i + 1] - p;
          r = (d[i] - g) * s + 2 * c * b;
          p = s * r;
          d[i + 1] = g + p;
          g = c * r - b;
          const ri = i * n, rj = (i + 1) * n;
          for (let k = 0; k < n; k++) {
            const f2 = zT[rj + k], z1 = zT[ri + k];
            zT[rj + k] = s * z1 + c * f2;
            zT[ri + k] = c * z1 - s * f2;
          }
        }
        if (r === 0 && i >= l) continue;
        d[l] -= p;
        e[l] = g;
        e[m] = 0;
      }
    } while (m !== l);
    if (onProgress && (l & 15) === 0) onProgress(l / n);
  }
}

export interface EigenResult {
  values: Float64Array; // ascending
  vectorsT: Float64Array; // (n, n): row k is the eigenvector of values[k]
}

/** All eigenpairs of the symmetric row-major matrix A (n x n, destroyed). */
export function symmetricEigen(A: Float64Array, n: number, onProgress?: (phase: string, frac: number) => void): EigenResult {
  const d = new Float64Array(n), e = new Float64Array(n);
  tred2(A, n, d, e, onProgress ? (f) => onProgress("tridiagonalisation", f) : undefined);
  // Q is in A (column j = j-th transformation); tql2 wants the transposed layout
  const zT = new Float64Array(n * n);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) zT[j * n + i] = A[i * n + j];
  tql2(d, e, n, zT, onProgress ? (f) => onProgress("QL iteration", f) : undefined);
  // sort ascending
  const order = Array.from({ length: n }, (_, i) => i).sort((p, q) => d[p] - d[q]);
  const values = new Float64Array(n), vectorsT = new Float64Array(n * n);
  order.forEach((k, i) => {
    values[i] = d[k];
    vectorsT.set(zT.subarray(k * n, (k + 1) * n), i * n);
  });
  return { values, vectorsT };
}