// bend-blas, the JS twins: plain loops, correct and slow, for the
// interpreter and the tests. An Array<F32> is a JS array of numbers.
// A failure is io_fail(code): the runtime's message for a code is not
// this library's, so the code is the contract on this lane.
function bend_fresh(n) {
  let slots = 1;
  while (slots < n) slots *= 2;
  return new Array(slots).fill(0);
}
function bend_copy(src, n) {
  const c = bend_fresh(n);
  for (let i = 0; i < n; i++) c[i] = src[i];
  return c;
}
function bend_short(a, n) {
  return n > 0x7fffffff || a.length < n;
}
function bend_gemm(ta, tb, m, n, k, alpha, a, b, beta, c) {
  const out = bend_copy(c, m * n);
  for (let i = 0; i < m; i++) for (let j = 0; j < n; j++) {
    let acc = 0;
    for (let l = 0; l < k; l++) acc += (ta ? a[l * m + i] : a[i * k + l]) * (tb ? b[j * k + l] : b[l * n + j]);
    out[i * n + j] = alpha * acc + beta * out[i * n + j];
  }
  return out;
}
function lapack_ssyev(n, a) {
  if (bend_short(a, n * n)) return io_fail(2);
  // cyclic Jacobi
  const A = bend_copy(a, n * n), V = new Array(n * n).fill(0);
  for (let i = 0; i < n; i++) { V[i * n + i] = 1; for (let j = i + 1; j < n; j++) A[i * n + j] = A[j * n + i]; }
  for (let sweep = 0; sweep < 60; sweep++) {
    let off = 0;
    for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) {
      const apq = A[p * n + q];
      if (Math.abs(apq) <= 1e-12) continue;
      off += Math.abs(apq);
      const theta = (A[q * n + q] - A[p * n + p]) / (2 * apq), t = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(1 + theta * theta));
      const cs = 1 / Math.sqrt(1 + t * t), sn = cs * t;
      for (let k = 0; k < n; k++) { const akp = A[k * n + p], akq = A[k * n + q]; A[k * n + p] = cs * akp - sn * akq; A[k * n + q] = sn * akp + cs * akq; }
      for (let k = 0; k < n; k++) { const apk = A[p * n + k], aqk = A[q * n + k]; A[p * n + k] = cs * apk - sn * aqk; A[q * n + k] = sn * apk + cs * aqk; }
      for (let k = 0; k < n; k++) { const vkp = V[k * n + p], vkq = V[k * n + q]; V[k * n + p] = cs * vkp - sn * vkq; V[k * n + q] = sn * vkp + cs * vkq; }
    }
    if (off === 0) break;
  }
  const order = [...Array(n).keys()].sort((x, y) => A[x * n + x] - A[y * n + y]);
  const w = bend_fresh(n), v = bend_fresh(n * n);
  for (let ii = 0; ii < n; ii++) { const i = order[ii]; w[ii] = A[i * n + i]; for (let k = 0; k < n; k++) v[ii * n + k] = V[k * n + i]; }
  return io_done({$: CID(Tuple), fst: w, snd: v});
}
io_eff(CID(Lapack.ssyev), lapack_ssyev);