// Eigenvalues of a symmetric tridiagonal matrix (diagonal d, off-diagonal e,
// e[i] couples i and i+1) by Sturm-sequence bisection. Zero off-diagonals
// split it into independent blocks, which the count handles.

function countBelow(d: Float64Array, e: Float64Array, x: number): number {
  let count = 0;
  let q = 1;
  for (let i = 0; i < d.length; i++) {
    const e2 = i > 0 ? e[i - 1] * e[i - 1] : 0;
    q = d[i] - x - (i > 0 ? e2 / q : 0);
    if (q === 0) q = -1e-300;
    if (q < 0) count++;
  }
  return count;
}

export function tridiagEigenvalues(d: Float64Array, e: Float64Array, count = d.length): number[] {
  const n = d.length;
  if (n === 0) return [];
  let lo = Infinity;
  let hi = -Infinity;
  for (let i = 0; i < n; i++) {
    const r = (i > 0 ? Math.abs(e[i - 1]) : 0) + (i < n - 1 ? Math.abs(e[i]) : 0);
    lo = Math.min(lo, d[i] - r);
    hi = Math.max(hi, d[i] + r);
  }
  const tol = 1e-13 * Math.max(Math.abs(lo), Math.abs(hi), 1e-300);
  const out: number[] = [];
  for (let k = 0; k < Math.min(count, n); k++) {
    let a = lo;
    let b = hi;
    while (b - a > tol) {
      const mid = 0.5 * (a + b);
      if (countBelow(d, e, mid) > k) b = mid;
      else a = mid;
      if (mid === a && mid === b) break;
    }
    out.push(0.5 * (a + b));
  }
  return out;
}