/**
 * MathEngine — Pure JS Linear Algebra Computation Library
 * Every method returns { value, steps: [{desc, latex}] }
 * Uses namespace pattern (not ES modules) for CDN compatibility.
 */
const MathEngine = (() => {
  const EPS = 0.001;

  // ── Helpers ──────────────────────────────────────────
  function round(x) {
    return Math.abs(x) < EPS ? 0 : Math.round(x * 10000) / 10000;
  }

  function cloneMatrix(M) {
    return M.map(row => [...row]);
  }

  function isMatrix(M) {
    return Array.isArray(M) && Array.isArray(M[0]);
  }

  function rows(M) { return M.length; }
  function cols(M) { return M[0].length; }

  function makeResult(value, steps) {
    return { value, steps: steps || [] };
  }

  function makeError(msg, steps) {
    return { error: msg, steps: steps || [] };
  }

  function vecToLatex(v) {
    return '(' + v.map(x => fmtNum(x)).join(';\\;') + ')';
  }

  function matToLatex(M) {
    const r = M.length;
    const c = M[0].length;
    let s = '\\begin{pmatrix}';
    for (let i = 0; i < r; i++) {
      s += M[i].map(x => fmtNum(x)).join(' & ');
      if (i < r - 1) s += ' \\\\ ';
    }
    s += '\\end{pmatrix}';
    return s;
  }

  function fmtNum(n) {
    if (n === undefined || n === null) return '0';
    const r = round(n);
    if (Number.isInteger(r)) return String(r);
    return r.toFixed(4).replace(/0+$/, '').replace(/\.$/, '');
  }

  // ── VECTOR OPERATIONS ───────────────────────────────
  const vector = {
    components(A, B) {
      const steps = [];
      const v = A.map((a, i) => round(B[i] - a));
      steps.push({ desc: 'Componentes del vector AB = B - A', latex: `\\overrightarrow{AB} = ${vecToLatex(B)} - ${vecToLatex(A)} = ${vecToLatex(v)}` });
      return makeResult(v, steps);
    },

    magnitude(v) {
      const steps = [];
      const squares = v.map(x => `${fmtNum(x)}^2`).join(' + ');
      const sum = v.reduce((s, x) => s + x * x, 0);
      const mag = round(Math.sqrt(sum));
      steps.push({ desc: 'M\u00f3dulo del vector', latex: `|\\vec{v}| = \\sqrt{${squares}} = \\sqrt{${fmtNum(sum)}} = ${fmtNum(mag)}` });
      return makeResult(mag, steps);
    },

    unitVector(v) {
      const steps = [];
      const mag = Math.sqrt(v.reduce((s, x) => s + x * x, 0));
      if (mag < EPS) return makeError('No se puede obtener vector unitario del vector nulo');
      const u = v.map(x => round(x / mag));
      steps.push({ desc: 'Vector unitario = v / |v|', latex: `\\hat{u} = \\frac{${vecToLatex(v)}}{${fmtNum(mag)}} = ${vecToLatex(u)}` });
      return makeResult(u, steps);
    },

    add(a, b) {
      const steps = [];
      const r = a.map((x, i) => round(x + b[i]));
      steps.push({ desc: 'Suma de vectores componente a componente', latex: `${vecToLatex(a)} + ${vecToLatex(b)} = ${vecToLatex(r)}` });
      return makeResult(r, steps);
    },

    scale(k, v) {
      const steps = [];
      const r = v.map(x => round(k * x));
      steps.push({ desc: 'Producto escalar', latex: `${fmtNum(k)} \\cdot ${vecToLatex(v)} = ${vecToLatex(r)}` });
      return makeResult(r, steps);
    },

    dotProduct(a, b) {
      const steps = [];
      const products = a.map((x, i) => `(${fmtNum(x)})(${fmtNum(b[i])})`).join(' + ');
      const val = a.reduce((s, x, i) => s + x * b[i], 0);
      const rv = round(val);
      steps.push({ desc: 'Producto escalar', latex: `\\vec{u} \\cdot \\vec{v} = ${products} = ${fmtNum(rv)}` });
      return makeResult(rv, steps);
    },

    crossProduct(a, b) {
      const steps = [];
      const r = [
        round(a[1] * b[2] - a[2] * b[1]),
        round(a[2] * b[0] - a[0] * b[2]),
        round(a[0] * b[1] - a[1] * b[0])
      ];
      steps.push({ desc: 'Producto vectorial por determinante', latex: `\\vec{u} \\times \\vec{v} = \\begin{vmatrix} \\mathbf{i} & \\mathbf{j} & \\mathbf{k} \\\\ ${fmtNum(a[0])} & ${fmtNum(a[1])} & ${fmtNum(a[2])} \\\\ ${fmtNum(b[0])} & ${fmtNum(b[1])} & ${fmtNum(b[2])} \\end{vmatrix}` });
      steps.push({ desc: 'Resultado', latex: `\\vec{u} \\times \\vec{v} = ${vecToLatex(r)}` });
      return makeResult(r, steps);
    },

    isParallel(a, b) {
      const steps = [];
      const cross = [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0]
      ];
      const isPar = cross.every(x => Math.abs(x) < EPS);
      steps.push({ desc: 'Son paralelos si u \u00d7 v = 0', latex: `\\vec{u} \\times \\vec{v} = ${vecToLatex(cross.map(round))} ${isPar ? '=' : '\\neq'} \\vec{0}` });
      steps.push({ desc: isPar ? 'Son paralelos' : 'No son paralelos', latex: `\\text{${isPar ? 'Son paralelos' : 'No son paralelos'}}` });
      return makeResult(isPar, steps);
    },

    isPerpendicular(a, b) {
      const steps = [];
      const dot = a.reduce((s, x, i) => s + x * b[i], 0);
      const isPerp = Math.abs(dot) < EPS;
      steps.push({ desc: 'Son perpendiculares si u \u00b7 v = 0', latex: `\\vec{u} \\cdot \\vec{v} = ${fmtNum(round(dot))} ${isPerp ? '=' : '\\neq'} 0` });
      steps.push({ desc: isPerp ? 'Son perpendiculares' : 'No son perpendiculares', latex: `\\text{${isPerp ? 'Son perpendiculares' : 'No son perpendiculares'}}` });
      return makeResult(isPerp, steps);
    },

    mixedProduct(a, b, c) {
      const steps = [];
      const cross = [
        b[1] * c[2] - b[2] * c[1],
        b[2] * c[0] - b[0] * c[2],
        b[0] * c[1] - b[1] * c[0]
      ];
      const val = a.reduce((s, x, i) => s + x * cross[i], 0);
      const rv = round(val);
      steps.push({ desc: 'Producto mixto [u,v,w] = u \u00b7 (v \u00d7 w)', latex: `\\vec{v} \\times \\vec{w} = ${vecToLatex(cross.map(round))}` });
      steps.push({ desc: 'Producto escalar con u', latex: `[\\vec{u},\\vec{v},\\vec{w}] = ${vecToLatex(a)} \\cdot ${vecToLatex(cross.map(round))} = ${fmtNum(rv)}` });
      return makeResult(rv, steps);
    },

    isCoplanar(a, b, c) {
      const steps = [];
      const mp = vector.mixedProduct(a, b, c);
      const isCo = Math.abs(mp.value) < EPS;
      steps.push(...mp.steps);
      steps.push({ desc: isCo ? 'Son coplanarios (producto mixto = 0)' : 'No son coplanarios (producto mixto \u2260 0)', latex: `[\\vec{u},\\vec{v},\\vec{w}] ${isCo ? '=' : '\\neq'} 0 \\Rightarrow \\text{${isCo ? 'Coplanarios' : 'No coplanarios'}}` });
      return makeResult(isCo, steps);
    },

    angle(a, b) {
      const steps = [];
      const dot = a.reduce((s, x, i) => s + x * b[i], 0);
      const magA = Math.sqrt(a.reduce((s, x) => s + x * x, 0));
      const magB = Math.sqrt(b.reduce((s, x) => s + x * x, 0));
      if (magA < EPS || magB < EPS) return makeError('Vector nulo no tiene \u00e1ngulo definido');
      const cosA = dot / (magA * magB);
      const clamped = Math.max(-1, Math.min(1, cosA));
      const angleRad = Math.acos(clamped);
      const angleDeg = round(angleRad * 180 / Math.PI);
      steps.push({ desc: '\u00c1ngulo entre vectores', latex: `\\cos\\alpha = \\frac{\\vec{u} \\cdot \\vec{v}}{|\\vec{u}||\\vec{v}|} = \\frac{${fmtNum(round(dot))}}{${fmtNum(round(magA))} \\cdot ${fmtNum(round(magB))}} = ${fmtNum(round(cosA))}` });
      steps.push({ desc: '\u00c1ngulo en grados', latex: `\\alpha = \\arccos(${fmtNum(round(cosA))}) = ${fmtNum(angleDeg)}^\\circ` });
      return makeResult(angleDeg, steps);
    }
  };

  // ── MATRIX OPERATIONS ───────────────────────────────
  const matrix = {
    create(rows, cols, fill = 0) {
      return Array.from({ length: rows }, () => Array(cols).fill(fill));
    },

    identity(n) {
      const M = matrix.create(n, n);
      for (let i = 0; i < n; i++) M[i][i] = 1;
      return M;
    },

    zeros(r, c) {
      if (c === undefined) c = r;
      return matrix.create(r, c, 0);
    },

    add(A, B) {
      const steps = [];
      const r = rows(A), c = cols(A);
      const R = matrix.create(r, c);
      for (let i = 0; i < r; i++)
        for (let j = 0; j < c; j++)
          R[i][j] = round(A[i][j] + B[i][j]);
      steps.push({ desc: 'Suma componente a componente', latex: `${matToLatex(A)} + ${matToLatex(B)} = ${matToLatex(R)}` });
      return makeResult(R, steps);
    },

    subtract(A, B) {
      const steps = [];
      const r = rows(A), c = cols(A);
      const R = matrix.create(r, c);
      for (let i = 0; i < r; i++)
        for (let j = 0; j < c; j++)
          R[i][j] = round(A[i][j] - B[i][j]);
      steps.push({ desc: 'Resta componente a componente', latex: `${matToLatex(A)} - ${matToLatex(B)} = ${matToLatex(R)}` });
      return makeResult(R, steps);
    },

    scale(k, A) {
      const steps = [];
      const R = A.map(row => row.map(x => round(k * x)));
      steps.push({ desc: 'Producto escalar por la matriz', latex: `${fmtNum(k)} \\cdot ${matToLatex(A)} = ${matToLatex(R)}` });
      return makeResult(R, steps);
    },

    multiply(A, B) {
      const steps = [];
      const rA = rows(A), cA = cols(A), cB = cols(B);
      if (cA !== B.length) return makeError(`Dimensiones incompatibles: ${rA}\u00d7${cA} por ${B.length}\u00d7${cB}`);
      const R = matrix.create(rA, cB);
      let detail = '';
      for (let i = 0; i < rA; i++) {
        for (let j = 0; j < cB; j++) {
          let sum = 0;
          const terms = [];
          for (let k = 0; k < cA; k++) {
            sum += A[i][k] * B[k][j];
            terms.push(`(${fmtNum(A[i][k])})(${fmtNum(B[k][j])})`);
          }
          R[i][j] = round(sum);
        }
      }
      steps.push({ desc: `Producto de matrices ${rA}\u00d7${cA} \u00b7 ${B.length}\u00d7${cB}`, latex: `${matToLatex(A)} \\cdot ${matToLatex(B)} = ${matToLatex(R)}` });
      return makeResult(R, steps);
    },

    transpose(A) {
      const steps = [];
      const r = rows(A), c = cols(A);
      const R = matrix.create(c, r);
      for (let i = 0; i < r; i++)
        for (let j = 0; j < c; j++)
          R[j][i] = A[i][j];
      steps.push({ desc: 'Traspuesta: intercambiar filas por columnas', latex: `(${matToLatex(A)})^T = ${matToLatex(R)}` });
      return makeResult(R, steps);
    },

    isSymmetric(A) {
      const n = rows(A);
      for (let i = 0; i < n; i++)
        for (let j = 0; j < i; j++)
          if (Math.abs(A[i][j] - A[j][i]) > EPS) return makeResult(false, [{ desc: 'No es sim\u00e9trica', latex: `a_{${i}${j}} = ${fmtNum(A[i][j])} \\neq a_{${j}${i}} = ${fmtNum(A[j][i])}` }]);
      return makeResult(true, [{ desc: 'La matriz es sim\u00e9trica', latex: 'A = A^T' }]);
    },

    trace(A) {
      const steps = [];
      const n = rows(A);
      const diag = [];
      let sum = 0;
      for (let i = 0; i < n; i++) {
        sum += A[i][i];
        diag.push(`a_{${i + 1}${i + 1}} = ${fmtNum(A[i][i])}`);
      }
      steps.push({ desc: 'Traza = suma de la diagonal', latex: `\\text{tr}(A) = ${diag.join(' + ')} = ${fmtNum(round(sum))}` });
      return makeResult(round(sum), steps);
    }
  };

  // ── DETERMINANT OPERATIONS ──────────────────────────
  const determinant = {
    det2x2(A) {
      const steps = [];
      const val = A[0][0] * A[1][1] - A[0][1] * A[1][0];
      const rv = round(val);
      steps.push({ desc: 'Determinante 2\u00d72', latex: `\\det ${matToLatex(A)} = (${fmtNum(A[0][0])})(${fmtNum(A[1][1])}) - (${fmtNum(A[0][1])})(${fmtNum(A[1][0])})` });
      steps.push({ desc: 'Resultado', latex: `= ${fmtNum(A[0][0] * A[1][1])} - ${fmtNum(A[0][1] * A[1][0])} = ${fmtNum(rv)}` });
      return makeResult(rv, steps);
    },

    det3x3Sarrus(A) {
      const steps = [];
      steps.push({ desc: 'Regla de Sarrus', latex: `\\det ${matToLatex(A)}` });
      const pos =
        A[0][0] * A[1][1] * A[2][2] +
        A[0][1] * A[1][2] * A[2][0] +
        A[0][2] * A[1][0] * A[2][1];
      const neg =
        A[0][2] * A[1][1] * A[2][0] +
        A[0][0] * A[1][2] * A[2][1] +
        A[0][1] * A[1][0] * A[2][2];
      const val = pos - neg;
      const rv = round(val);
      steps.push({ desc: 'Diagonales positivas', latex: `+${fmtNum(A[0][0])} \\cdot ${fmtNum(A[1][1])} \\cdot ${fmtNum(A[2][2])} + ${fmtNum(A[0][1])} \\cdot ${fmtNum(A[1][2])} \\cdot ${fmtNum(A[2][0])} + ${fmtNum(A[0][2])} \\cdot ${fmtNum(A[1][0])} \\cdot ${fmtNum(A[2][1])}` });
      steps.push({ desc: 'Suma positivas', latex: `= ${fmtNum(round(pos))}` });
      steps.push({ desc: 'Diagonales negativas', latex: `-${fmtNum(A[0][2])} \\cdot ${fmtNum(A[1][1])} \\cdot ${fmtNum(A[2][0])} - ${fmtNum(A[0][0])} \\cdot ${fmtNum(A[1][2])} \\cdot ${fmtNum(A[2][1])} - ${fmtNum(A[0][1])} \\cdot ${fmtNum(A[1][0])} \\cdot ${fmtNum(A[2][2])}` });
      steps.push({ desc: 'Suma negativas', latex: `= -${fmtNum(round(neg))}` });
      steps.push({ desc: 'Resultado', latex: `${fmtNum(round(pos))} - ${fmtNum(round(neg))} = ${fmtNum(rv)}` });
      return makeResult(rv, steps);
    },

    cofactor(A, row, col) {
      const n = rows(A);
      const sub = [];
      for (let i = 0; i < n; i++) {
        if (i === row) continue;
        const r = [];
        for (let j = 0; j < n; j++) {
          if (j === col) continue;
          r.push(A[i][j]);
        }
        sub.push(r);
      }
      const sign = (row + col) % 2 === 0 ? 1 : -1;
      const subDet = determinant.det(sub);
      const val = sign * subDet.value;
      const rv = round(val);
      const fullSteps = [];
      fullSteps.push({ desc: `Cofactor C_{${row + 1}${col + 1}}`, latex: `C_{${row + 1}${col + 1}} = ${sign > 0 ? '' : '-'} \\det ${matToLatex(sub)}` });
      fullSteps.push(...subDet.steps);
      fullSteps.push({ desc: 'Valor del cofactor', latex: `C_{${row + 1}${col + 1}} = ${fmtNum(rv)}` });
      return makeResult(rv, fullSteps);
    },

    det(A) {
      const n = rows(A);
      if (n === 1) return makeResult(A[0][0], [{ desc: 'Determinante 1\u00d71', latex: `\\det = ${fmtNum(A[0][0])}` }]);
      if (n === 2) return determinant.det2x2(A);
      if (n === 3) return determinant.det3x3Sarrus(A);
      // Generic cofactor expansion along first row
      const steps = [];
      steps.push({ desc: `Expansi\u00f3n por cofactores (fila 1)`, latex: `\\det ${matToLatex(A)}` });
      let total = 0;
      const parts = [];
      for (let j = 0; j < n; j++) {
        const sign = j % 2 === 0 ? 1 : -1;
        const subDet = determinant.cofactor(A, 0, j);
        const contrib = sign * A[0][j] * subDet.value;
        total += contrib;
        parts.push(`${sign > 0 ? '' : '-'} ${fmtNum(A[0][j])} \\cdot C_{1${j + 1}}`);
        steps.push(...subDet.steps);
      }
      steps.push({ desc: 'Expansi\u00f3n completa', latex: parts.join(' + ') });
      steps.push({ desc: 'Resultado', latex: `= ${fmtNum(round(total))}` });
      return makeResult(round(total), steps);
    },

    detByTriangularization(A) {
      const steps = [];
      const n = rows(A);
      const M = cloneMatrix(A);
      let sign = 1;
      steps.push({ desc: 'Triangularizaci\u00f3n', latex: `\\det ${matToLatex(M)}` });
      for (let col = 0; col < n; col++) {
        // Find pivot
        let pivotRow = -1;
        for (let i = col; i < n; i++) {
          if (Math.abs(M[i][col]) > EPS) { pivotRow = i; break; }
        }
        if (pivotRow === -1) {
          steps.push({ desc: 'Columna sin pivote \u2192 determinante 0', latex: `\\det = 0` });
          return makeResult(0, steps);
        }
        if (pivotRow !== col) {
          [M[col], M[pivotRow]] = [M[pivotRow], M[col]];
          sign *= -1;
          steps.push({ desc: `F${col + 1} \u2194 F${pivotRow + 1} (cambia signo)`, latex: matToLatex(M) });
        }
        for (let i = col + 1; i < n; i++) {
          if (Math.abs(M[i][col]) > EPS) {
            const factor = M[i][col] / M[col][col];
            for (let j = col; j < n; j++) {
              M[i][j] = round(M[i][j] - factor * M[col][j]);
            }
            steps.push({ desc: `F${i + 1} \u2192 F${i + 1} - (${fmtNum(round(factor))})F${col + 1}`, latex: matToLatex(M) });
          }
        }
      }
      let detVal = sign;
      for (let i = 0; i < n; i++) detVal *= M[i][i];
      detVal = round(detVal);
      const diagParts = [];
      for (let i = 0; i < n; i++) diagParts.push(fmtNum(M[i][i]));
      steps.push({ desc: 'Producto de la diagonal', latex: `\\det = ${sign < 0 ? '-' : ''}${diagParts.join(' \\cdot ')} = ${fmtNum(detVal)}` });
      return makeResult(detVal, steps);
    }
  };

  // ── INVERSE OPERATIONS ──────────────────────────────
  const inverse = {
    inverse2x2(A) {
      const steps = [];
      const detResult = determinant.det2x2(A);
      steps.push(...detResult.steps);
      const det = detResult.value;
      if (Math.abs(det) < EPS) {
        steps.push({ desc: 'Matriz no invertible (det = 0)', latex: '\\det = 0 \\Rightarrow A^{-1} \\text{ no existe}' });
        return makeError('Matriz no invertible (determinante nulo)', steps);
      }
      const R = [[A[1][1] / det, -A[0][1] / det], [-A[1][0] / det, A[0][0] / det]].map(r => r.map(x => round(x)));
      steps.push({ desc: 'Inversa 2\u00d72: intercambiar diagonal, negar antidiagonal, dividir por det', latex: `A^{-1} = \\frac{1}{${fmtNum(det)}} \\begin{pmatrix} ${fmtNum(A[1][1])} & ${fmtNum(-A[0][1])} \\\\ ${fmtNum(-A[1][0])} & ${fmtNum(A[0][0])} \\end{pmatrix}` });
      steps.push({ desc: 'Resultado', latex: `A^{-1} = ${matToLatex(R)}` });
      return makeResult(R, steps);
    },

    inverse(A) {
      const n = rows(A);
      if (n === 2) return inverse.inverse2x2(A);
      const steps = [];
      const detResult = determinant.det(A);
      steps.push(...detResult.steps);
      const det = detResult.value;
      if (Math.abs(det) < EPS) {
        steps.push({ desc: 'Matriz no invertible (det = 0)', latex: '\\det = 0 \\Rightarrow A^{-1} \\text{ no existe}' });
        return makeError('Matriz no invertible (determinante nulo)', steps);
      }
      // Build adjugate (transpose of cofactor matrix)
      steps.push({ desc: 'Construir la adjunta (traspuesta de la matriz de cofactores)', latex: `\\text{Adj}(A)` });
      const adj = matrix.create(n, n);
      for (let i = 0; i < n; i++) {
        for (let j = 0; j < n; j++) {
          const cf = determinant.cofactor(A, i, j);
          adj[j][i] = cf.value; // note: transposing
        }
      }
      const R = adj.map(row => row.map(x => round(x / det)));
      steps.push({ desc: 'Dividir adjunta por determinante', latex: `A^{-1} = \\frac{1}{${fmtNum(det)}} \\cdot ${matToLatex(adj)}` });
      steps.push({ desc: 'Resultado', latex: `A^{-1} = ${matToLatex(R)}` });
      return makeResult(R, steps);
    },

    isInvertible(A) {
      const det = determinant.det(A);
      const inv = Math.abs(det.value) > EPS;
      return makeResult(inv, [...det.steps, { desc: inv ? 'Es invertible (det \u2260 0)' : 'No es invertible (det = 0)', latex: inv ? '\\det(A) \\neq 0 \\Rightarrow \\text{Invertible}' : '\\det(A) = 0 \\Rightarrow \\text{No invertible}' }]);
    }
  };

  // ── SYSTEM SOLVERS ──────────────────────────────────
  const system = {
    rank(A) {
      const steps = [];
      const M = cloneMatrix(A);
      const r = rows(M), c = cols(M);
      const n = Math.min(r, c);
      let rk = 0;
      steps.push({ desc: 'Reducir a forma escalonada', latex: matToLatex(M) });
      for (let col = 0; col < c && rk < r; col++) {
        let pivotRow = -1;
        for (let i = rk; i < r; i++) {
          if (Math.abs(M[i][col]) > EPS) { pivotRow = i; break; }
        }
        if (pivotRow === -1) continue;
        if (pivotRow !== rk) {
          [M[rk], M[pivotRow]] = [M[pivotRow], M[rk]];
          steps.push({ desc: `F${rk + 1} \u2194 F${pivotRow + 1}`, latex: matToLatex(M) });
        }
        const pivot = M[rk][col];
        for (let i = rk + 1; i < r; i++) {
          if (Math.abs(M[i][col]) > EPS) {
            const factor = M[i][col] / pivot;
            for (let j = col; j < c; j++) {
              M[i][j] = round(M[i][j] - factor * M[rk][j]);
            }
            steps.push({ desc: `F${i + 1} \u2192 F${i + 1} - (${fmtNum(round(factor))})F${rk + 1}`, latex: matToLatex(M) });
          }
        }
        rk++;
      }
      steps.push({ desc: `Rango = ${rk}`, latex: `\\text{rg}(A) = ${rk}` });
      return makeResult(rk, steps);
    },

    gaussElimination(A, b) {
      const steps = [];
      const n = A.length;
      const aug = A.map((row, i) => [...row, b[i]]);

      steps.push({ desc: 'Sistema aumentado inicial', latex: matToLatex(aug) });

      // Forward elimination
      for (let col = 0; col < n; col++) {
        let pivotRow = -1;
        for (let i = col; i < n; i++) {
          if (Math.abs(aug[i][col]) > EPS) { pivotRow = i; break; }
        }
        if (pivotRow === -1) {
          // Check for incompatibility
          for (let i = col; i < n; i++) {
            if (Math.abs(aug[i][n]) > EPS) {
              steps.push({ desc: 'Sistema incompatible', latex: `0 = ${fmtNum(round(aug[i][n]))} \\Rightarrow \\text{Incompatible}` });
              return makeError('Sistema incompatible', steps);
            }
          }
          continue;
        }
        if (pivotRow !== col) {
          [aug[col], aug[pivotRow]] = [aug[pivotRow], aug[col]];
          steps.push({ desc: `F${col + 1} \u2194 F${pivotRow + 1}`, latex: matToLatex(aug) });
        }
        for (let i = col + 1; i < n; i++) {
          if (Math.abs(aug[i][col]) > EPS) {
            const factor = aug[i][col] / aug[col][col];
            for (let j = col; j <= n; j++) {
              aug[i][j] = round(aug[i][j] - factor * aug[col][j]);
            }
            steps.push({ desc: `F${i + 1} \u2192 F${i + 1} - (${fmtNum(round(factor))})F${col + 1}`, latex: matToLatex(aug) });
          }
        }
      }

      // Check for indeterminate
      const rankA = aug.filter(row => row.slice(0, n).some(x => Math.abs(x) > EPS)).length;
      if (rankA < n) {
        steps.push({ desc: `Sistema compatible indeterminado (rango = ${rankA} < ${n})`, latex: `\\text{rg}(A) = ${rankA} < ${n} \\Rightarrow \\text{CI}` });
        return makeResult({ type: 'indeterminate', rank: rankA, augmented: aug }, steps);
      }

      // Back substitution
      const x = new Array(n);
      for (let i = n - 1; i >= 0; i--) {
        let sum = aug[i][n];
        for (let j = i + 1; j < n; j++) {
          sum -= aug[i][j] * x[j];
        }
        x[i] = round(sum / aug[i][i]);
      }

      steps.push({ desc: 'Sustituci\u00f3n regresiva', latex: `x = ${vecToLatex(x)}` });
      return makeResult(x, steps);
    },

    gaussJordan(A, b) {
      const steps = [];
      const n = A.length;
      const aug = A.map((row, i) => [...row, b[i]]);
      steps.push({ desc: 'Matriz aumentada', latex: matToLatex(aug) });

      // Forward elimination
      let pivotCols = [];
      let rowIdx = 0;
      for (let col = 0; col < n && rowIdx < n; col++) {
        let pivotRow = -1;
        for (let i = rowIdx; i < n; i++) {
          if (Math.abs(aug[i][col]) > EPS) { pivotRow = i; break; }
        }
        if (pivotRow === -1) continue;
        if (pivotRow !== rowIdx) {
          [aug[rowIdx], aug[pivotRow]] = [aug[pivotRow], aug[rowIdx]];
          steps.push({ desc: `F${rowIdx + 1} \u2194 F${pivotRow + 1}`, latex: matToLatex(aug) });
        }
        const pivot = aug[rowIdx][col];
        for (let j = col; j <= n; j++) aug[rowIdx][j] = round(aug[rowIdx][j] / pivot);
        steps.push({ desc: `F${rowIdx + 1} \u2192 F${rowIdx + 1} / ${fmtNum(round(pivot))}`, latex: matToLatex(aug) });
        for (let i = 0; i < n; i++) {
          if (i !== rowIdx && Math.abs(aug[i][col]) > EPS) {
            const factor = aug[i][col];
            for (let j = col; j <= n; j++) {
              aug[i][j] = round(aug[i][j] - factor * aug[rowIdx][j]);
            }
            steps.push({ desc: `F${i + 1} \u2192 F${i + 1} - (${fmtNum(round(factor))})F${rowIdx + 1}`, latex: matToLatex(aug) });
          }
        }
        pivotCols.push(col);
        rowIdx++;
      }

      const rankA = pivotCols.length;
      // Check for incompatibility
      for (let i = rankA; i < n; i++) {
        if (Math.abs(aug[i][n]) > EPS) {
          steps.push({ desc: 'Sistema incompatible', latex: `0 = ${fmtNum(round(aug[i][n]))}` });
          return makeError('Sistema incompatible (SI)', steps);
        }
      }
      if (rankA < n) {
        steps.push({ desc: `Sistema compatible indeterminado (rango = ${rankA})`, latex: `\\text{rg}(A) = ${rankA} < ${n}` });
        const x = new Array(n).fill(0);
        for (let i = 0; i < rankA; i++) x[pivotCols[i]] = round(aug[i][n]);
        return makeResult({ type: 'indeterminate', rank: rankA, solution: x, augmented: aug, pivotCols }, steps);
      }

      const x = aug.map(row => round(row[n]));
      steps.push({ desc: 'Soluci\u00f3n', latex: `x = ${vecToLatex(x)}` });
      return makeResult(x, steps);
    },

    cramer(A, b) {
      const steps = [];
      const n = A.length;
      const detA = determinant.det(A);
      steps.push({ desc: 'Determinante del sistema', latex: `\\Delta = \\det(A) = ${fmtNum(detA.value)}` });
      steps.push(...detA.steps);

      if (Math.abs(detA.value) < EPS) {
        steps.push({ desc: 'No se puede aplicar Cramer (det = 0)', latex: '\\Delta = 0 \\Rightarrow \\text{Cramer no aplica}' });
        return makeError('Sistema con determinante nulo (Cramer no aplica)', steps);
      }

      const x = [];
      for (let j = 0; j < n; j++) {
        const Aj = cloneMatrix(A);
        for (let i = 0; i < n; i++) Aj[i][j] = b[i];
        const detJ = determinant.det(Aj);
        steps.push({ desc: `\\Delta_${j + 1}: reemplazar columna ${j + 1} por t\u00e9rminos independientes`, latex: `\\Delta_${j + 1} = \\det ${matToLatex(Aj)} = ${fmtNum(detJ.value)}` });
        x.push(round(detJ.value / detA.value));
      }

      const vars = x.map((v, i) => `x_{${i + 1}} = \\frac{${fmtNum(round(determinant.det(A.map((row, ii) => { const Aj = cloneMatrix(A); for (let k = 0; k < n; k++) Aj[k][i] = b[k]; return Aj; })).value))}}{${fmtNum(detA.value)}} = ${fmtNum(v)}`).join(', \\; ');
      steps.push({ desc: 'Soluci\u00f3n por Cramer', latex: vars });
      steps.push({ desc: 'Resultado', latex: `x = ${vecToLatex(x)}` });
      return makeResult(x, steps);
    },

    roucheFrobenius(A, b) {
      const steps = [];
      const n = A.length;
      const augMatrix = A.map((row, i) => [...row, b[i]]);
      const rankA = system.rank(cloneMatrix(A));
      const rankAug = system.rank(cloneMatrix(augMatrix));
      steps.push(...rankA.steps);
      steps.push(...rankAug.steps);

      const rA = rankA.value;
      const rAug = rankAug.value;

      let classification;
      if (rA !== rAug) {
        classification = { type: 'SI', label: 'Sistema Incompatible' };
      } else if (rA === rAug && rA === n) {
        classification = { type: 'CD', label: 'Compatible Determinado' };
      } else {
        classification = { type: 'CI', label: 'Compatible Indeterminado', freeVars: n - rA };
      }

      steps.push({ desc: 'Clasificaci\u00f3n Rouche-Frobenius', latex: `\\text{rg}(A) = ${rA}, \\; \\text{rg}(A|b) = ${rAug}, \\; n = ${n} \\Rightarrow ${classification.label}` });
      return makeResult(classification, steps);
    },

    solveHomogeneous(A) {
      const steps = [];
      const n = A.length;
      const b = new Array(n).fill(0);
      steps.push({ desc: 'Sistema homog\u00e9neo Ax = 0', latex: matToLatex(A.map((row, i) => [...row, 0])) });

      const rankA = system.rank(cloneMatrix(A));
      steps.push(...rankA.steps);

      if (rankA.value === n) {
        steps.push({ desc: 'Rango = n, soluci\u00f3n \u00fanica trivial', latex: `\\text{rg}(A) = ${n} = n \\Rightarrow x = \\vec{0}` });
        return makeResult({ type: 'trivial', solution: new Array(n).fill(0) }, steps);
      }

      const gj = system.gaussJordan(cloneMatrix(A), [...b]);
      steps.push(...gj.steps);

      if (gj.error) {
        return makeResult({ type: 'trivial', solution: new Array(n).fill(0) }, steps);
      }

      if (gj.value && gj.value.type === 'indeterminate') {
        steps.push({ desc: `Soluci\u00f3n no trivial (${n - rankA.value} par\u00e1metros libres)`, latex: `\\dim(\\ker) = ${n - rankA.value}` });
      }

      return makeResult(gj.value || { type: 'trivial', solution: new Array(n).fill(0) }, steps);
    }
  };

  return {
    EPS,
    round,
    cloneMatrix,
    vecToLatex,
    matToLatex,
    fmtNum,
    vector,
    matrix,
    determinant,
    inverse,
    system
  };
})();