// 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_sgesvd_s(m, n, a) {
  if (bend_short(a, m * n)) return io_fail(2);
  // one-sided Jacobi on the taller orientation
  const tall = m >= n, R = tall ? m : n, C = tall ? n : m;
  const W = new Array(R * C);
  for (let i = 0; i < R; i++) for (let j = 0; j < C; j++) W[i * C + j] = tall ? a[i * n + j] : a[j * n + i];
  const V = new Array(C * C).fill(0);
  for (let i = 0; i < C; i++) V[i * C + i] = 1;
  for (let sweep = 0; sweep < 60; sweep++) {
    let off = 0;
    for (let p = 0; p < C; p++) for (let q = p + 1; q < C; q++) {
      let al = 0, be = 0, ga = 0;
      for (let i = 0; i < R; i++) { al += W[i * C + p] * W[i * C + p]; be += W[i * C + q] * W[i * C + q]; ga += W[i * C + p] * W[i * C + q]; }
      if (Math.abs(ga) <= 1e-12 * Math.sqrt(al * be)) continue;
      off += Math.abs(ga);
      const zeta = (be - al) / (2 * ga), t = Math.sign(zeta || 1) / (Math.abs(zeta) + Math.sqrt(1 + zeta * zeta));
      const cs = 1 / Math.sqrt(1 + t * t), sn = cs * t;
      for (let i = 0; i < R; i++) { const wp = W[i * C + p], wq = W[i * C + q]; W[i * C + p] = cs * wp - sn * wq; W[i * C + q] = sn * wp + cs * wq; }
      for (let i = 0; i < C; i++) { const vp = V[i * C + p], vq = V[i * C + q]; V[i * C + p] = cs * vp - sn * vq; V[i * C + q] = sn * vp + cs * vq; }
    }
    if (off === 0) break;
  }
  const order = [...Array(C).keys()], norm = new Array(C);
  for (let j = 0; j < C; j++) { let s = 0; for (let i = 0; i < R; i++) s += W[i * C + j] * W[i * C + j]; norm[j] = Math.sqrt(s); }
  order.sort((x, y) => norm[y] - norm[x]);
  const k = C, s = bend_fresh(k), Ut = new Array(R * k), Vt = new Array(k * C);
  for (let jj = 0; jj < k; jj++) {
    const j = order[jj]; s[jj] = norm[j];
    for (let i = 0; i < R; i++) Ut[i * k + jj] = norm[j] > 0 ? W[i * C + j] / norm[j] : 0;
    for (let i = 0; i < C; i++) Vt[jj * C + i] = V[i * C + j];
  }
  // for a wide input the roles swap: A = (A^T)^T = V S U^T
  const u = bend_fresh(m * k), vt = bend_fresh(k * n);
  for (let i = 0; i < m; i++) for (let jj = 0; jj < k; jj++) u[i * k + jj] = tall ? Ut[i * k + jj] : Vt[jj * C + i];
  for (let jj = 0; jj < k; jj++) for (let j = 0; j < n; j++) vt[jj * n + j] = tall ? Vt[jj * C + j] : Ut[j * k + jj];
  return io_done(s);
}
io_eff(CID(Lapack.sgesvd_s), lapack_sgesvd_s);