///|
fn swap_dense_rows(data : Array[Array[Expr]], i : Int, j : Int) -> Unit {
  let tmp = data[i]
  data[i] = data[j]
  data[j] = tmp
}

///|
fn symbolic_pivot_better(candidate : Expr, current : Expr) -> Bool {
  let candidate_exact = @symcore.exact_number_num_den(candidate) is Some(_)
  let current_exact = @symcore.exact_number_num_den(current) is Some(_)
  if candidate_exact != current_exact {
    return candidate_exact
  }
  let candidate_free = @symcore.free_symbols(candidate).length()
  let current_free = @symcore.free_symbols(current).length()
  if candidate_free != current_free {
    return candidate_free < current_free
  }
  @symprint.pretty_string(candidate).length() <
  @symprint.pretty_string(current).length()
}

///|
fn symbolic_pivot_row(
  data : Array[Array[Expr]],
  row : Int,
  rows : Int,
  col : Int,
) -> Int {
  let mut pivot_row = -1
  for candidate in row.. Expr {
  @symcore.add([lhs, @symcore.mul([@symcore.int(-1), rhs])])
}

///|
fn raw_expr_mul(lhs : Expr, rhs : Expr) -> Expr {
  @symcore.mul([lhs, rhs])
}

///|
fn raw_expr_div(lhs : Expr, rhs : Expr) -> Expr {
  @symcore.mul([lhs, @symcore.pow(rhs, @symcore.int(-1))])
}

///|
fn matrix_exact_dm(
  rows : Array[Array[Expr]],
) -> @sympolys.DomainMatrix raise MatrixError {
  @sympolys.DomainMatrix::from_Matrix(rows, fmt="dense") catch {
    _ =>
      raise MatrixError::ValueError(
        "matrix entries are not supported by the exact DomainMatrix linear algebra path",
      )
  }
}

///|
fn matrix_from_domain_matrix(
  dm : @sympolys.DomainMatrix,
) -> Matrix raise MatrixError {
  let rows = dm.to_Matrix() catch {
    _ =>
      raise MatrixError::ValueError(
        "exact DomainMatrix result could not be converted back to a symbolic matrix",
      )
  }
  { rows: dm.rows, cols: dm.cols, data: rows }
}

///|
fn matrix_index_range(start : Int, end_ : Int) -> Array[Int] {
  let out : Array[Int] = Array::new()
  for i in start.. Matrix raise MatrixError {
  match matrix_exact_integer_data(matrix.to_list()) {
    Some(data) =>
      match exact_integer_inverse(data) {
        Some((num, den)) =>
          return {
            rows: matrix.rows,
            cols: matrix.cols,
            data: exact_integer_fraction_rows_to_expr_rows(num, den),
          }
        None => raise MatrixError::SingularMatrixError("matrix is singular")
      }
    None => ()
  }
  let dm = matrix_exact_dm(matrix.to_list())
  let identity = @sympolys.DomainMatrix::eye(matrix.rows, dm.domain) catch {
    _ =>
      raise MatrixError::ValueError(
        "exact DomainMatrix inverse identity construction failed",
      )
  }
  let aug = dm.hstack(identity) catch {
    _ =>
      raise MatrixError::ValueError(
        "exact DomainMatrix inverse augmentation failed",
      )
  }
  let (aug_rref, den0, pivots) = aug.rref_den(rref_method="FF_dense") catch {
    _ => raise MatrixError::SingularMatrixError("matrix is singular")
  }
  if pivots.length() != matrix.rows {
    raise MatrixError::SingularMatrixError("matrix is singular")
  }
  for i in 0..
      raise MatrixError::ValueError(
        "exact DomainMatrix inverse extraction failed",
      )
  }
  let num = if @sympolys.domain_is_field_pred(num0.domain) {
    num0
  } else {
    num0.to_field() catch {
      _ =>
        raise MatrixError::ValueError(
          "exact DomainMatrix inverse could not be converted into its field domain",
        )
    }
  }
  let den = (@sympolys.DomainScalar::new(den0, num0.domain) catch {
      _ =>
        raise MatrixError::ValueError(
          "exact DomainMatrix inverse denominator could not be reconstructed",
        )
    }).convert_to(num.domain).element catch {
      _ =>
        raise MatrixError::ValueError(
          "exact DomainMatrix inverse denominator could not be converted into the field domain",
        )
    }
  let inverse = num.mul(@sympolys.fe_inv(den)) catch {
    _ =>
      raise MatrixError::ValueError(
        "exact DomainMatrix inverse could not be normalized in its field domain",
      )
  }
  matrix_from_domain_matrix(inverse)
}

///|
fn matrix_exact_solve(lhs : Matrix, rhs : Matrix) -> Matrix raise MatrixError {
  match
    (
      matrix_exact_integer_data(lhs.to_list()),
      matrix_exact_integer_data(rhs.to_list()),
    ) {
    (Some(lhs_data), Some(rhs_data)) =>
      match exact_integer_solve(lhs_data, rhs_data) {
        Some((num, den)) =>
          return {
            rows: lhs.cols,
            cols: rhs.cols,
            data: exact_integer_fraction_rows_to_expr_rows(num, den),
          }
        None => raise MatrixError::SingularMatrixError("matrix is singular")
      }
    _ => ()
  }
  let lhs_dm = matrix_exact_dm(lhs.to_list())
  let rhs_dm = matrix_exact_dm(rhs.to_list())
  let aug = lhs_dm.hstack(rhs_dm) catch {
    _ =>
      raise MatrixError::ValueError(
        "exact DomainMatrix solve augmentation failed",
      )
  }
  let (aug_rref, den0, pivots) = aug.rref_den(rref_method="FF_dense") catch {
    _ => raise MatrixError::SingularMatrixError("matrix is singular")
  }
  if pivots.length() != lhs.cols {
    raise MatrixError::SingularMatrixError("matrix is singular")
  }
  if !pivots.is_empty() && pivots[pivots.length() - 1] >= lhs.cols {
    raise MatrixError::SingularMatrixError("matrix is singular")
  }
  let num0 = aug_rref.extract(
    matrix_index_range(0, lhs.cols),
    matrix_index_range(lhs.cols, aug.cols),
  ) catch {
    _ =>
      raise MatrixError::ValueError(
        "exact DomainMatrix solve extraction failed",
      )
  }
  let num = if @sympolys.domain_is_field_pred(num0.domain) {
    num0
  } else {
    num0.to_field() catch {
      _ =>
        raise MatrixError::ValueError(
          "exact DomainMatrix solve result could not be converted into its field domain",
        )
    }
  }
  let den = (@sympolys.DomainScalar::new(den0, num0.domain) catch {
      _ =>
        raise MatrixError::ValueError(
          "exact DomainMatrix solve denominator could not be reconstructed",
        )
    }).convert_to(num.domain).element catch {
      _ =>
        raise MatrixError::ValueError(
          "exact DomainMatrix solve denominator could not be converted into the field domain",
        )
    }
  let solution = num.mul(@sympolys.fe_inv(den)) catch {
    _ =>
      raise MatrixError::ValueError(
        "exact DomainMatrix solve result could not be normalized in its field domain",
      )
  }
  matrix_from_domain_matrix(solution)
}

///|
fn matrix_try_exact_inverse(matrix : Matrix) -> Matrix? raise MatrixError {
  let exact : Result[Matrix, MatrixError] = try? matrix_exact_inverse(matrix)
  match exact {
    Ok(value) => Some(value)
    Err(MatrixError::ValueError(_)) => None
    Err(err) => raise err
  }
}

///|
fn matrix_try_exact_solve(
  lhs : Matrix,
  rhs : Matrix,
) -> Matrix? raise MatrixError {
  let exact : Result[Matrix, MatrixError] = try? matrix_exact_solve(lhs, rhs)
  match exact {
    Ok(value) => Some(value)
    Err(MatrixError::ValueError(_)) => None
    Err(err) => raise err
  }
}

///|
fn matrix_echelon_internal(matrix : Matrix) -> (Array[Array[Expr]], Array[Int]) {
  let data = clone_expr_rows(matrix.data)
  let pivots : Array[Int] = []
  let mut row = 0
  let mut col = 0
  let mut prev = @symcore.int(1)
  while row < matrix.rows && col < matrix.cols {
    let pivot_row = symbolic_pivot_row(data, row, matrix.rows, col)
    if pivot_row == -1 {
      col += 1
      continue
    }
    if pivot_row != row {
      swap_dense_rows(data, pivot_row, row)
    }
    let pivot = data[row][col]
    for i in (row + 1).. (Array[Array[Expr]], Array[Int]) {
  let data = clone_expr_rows(matrix.data)
  let pivots : Array[Int] = []
  let mut row = 0
  let mut col = 0
  while row < matrix.rows && col < matrix.cols {
    let pivot_row = symbolic_pivot_row(data, row, matrix.rows, col)
    if pivot_row == -1 {
      col += 1
      continue
    }
    if pivot_row != row {
      swap_dense_rows(data, pivot_row, row)
    }
    let pivot = data[row][col]
    for j in col.. (Matrix, Array[Int]) {
  let numeric_prec = expr_numeric_precision(self.data)
  if matrix_has_floating_leaf(self.data) &&
    matrix_all_numeric(self.data, numeric_prec) {
    let data = clone_expr_rows(self.data).map(row => {
      row.map(cell => @symcore.evalf(cell, prec=numeric_prec))
    })
    let pivots : Array[Int] = []
    let mut row = 0
    let mut col = 0
    while row < self.rows && col < self.cols {
      let mut pivot_row = row
      let mut pivot_measure = expr_numeric_abs(data[row][col], numeric_prec)
      for candidate in (row + 1).. 0 => {
            pivot_row = candidate
            pivot_measure = Some(measure)
          }
          (Some(measure), None) => {
            pivot_row = candidate
            pivot_measure = Some(measure)
          }
          _ => ()
        }
      }
      if expr_numeric_zero(data[pivot_row][col], numeric_prec) {
        col += 1
        continue
      }
      if pivot_row != row {
        swap_dense_rows(data, pivot_row, row)
      }
      let pivot = data[row][col]
      for j in col.. Int {
  self.rref().1.length()
}

///|
/// Compute an echelon form of a matrix.
///
/// - Does: Performs elimination without normalizing pivots to one.
/// - Input: Any dense `Matrix`.
/// - Returns: A dense `Matrix` in echelon form.
/// - Limits: This front door does not report pivot columns; use `echelon_form_with_pivots` when callers need them.
pub fn Matrix::echelon_form(self : Matrix) -> Matrix {
  let (data, _) = matrix_echelon_internal(self)
  { rows: self.rows, cols: self.cols, data }
}

///|
pub fn Matrix::echelon_form_with_pivots(self : Matrix) -> (Matrix, Array[Int]) {
  let (data, pivots) = matrix_echelon_internal(self)
  ({ rows: self.rows, cols: self.cols, data }, pivots)
}

///|
pub fn Matrix::is_echelon(self : Matrix) -> Bool {
  let mut last_pivot = -1
  let mut seen_zero_row = false
  for i in 0.. Expr raise MatrixError {
  if !self.is_square() {
    raise MatrixError::NonSquareMatrixError(
      "determinant requires a square matrix",
    )
  }
  if self.rows == 0 {
    return @symcore.int(1)
  }
  match matrix_exact_rational_data(self.data) {
    Some(data) => return exact_expr_from_rational(exact_rational_det(data))
    None => ()
  }
  let numeric_prec = expr_numeric_precision(self.data)
  if matrix_has_floating_leaf(self.data) &&
    matrix_all_numeric(self.data, numeric_prec) {
    let data = clone_expr_rows(self.evalf(prec=numeric_prec).data)
    let mut sign = @symcore.int(1)
    let mut col = 0
    while col < self.rows {
      let mut pivot_row = col
      let mut pivot_measure = expr_numeric_abs(data[col][col], numeric_prec)
      for candidate in (col + 1).. 0 => {
            pivot_row = candidate
            pivot_measure = Some(measure)
          }
          (Some(measure), None) => {
            pivot_row = candidate
            pivot_measure = Some(measure)
          }
          _ => ()
        }
      }
      if expr_numeric_zero(data[pivot_row][col], numeric_prec) {
        return @symcore.int(0)
      }
      if pivot_row != col {
        swap_dense_rows(data, pivot_row, col)
        sign = expr_numeric_mul(sign, expr_neg(@symcore.int(1)), numeric_prec)
      }
      let pivot = data[col][col]
      for row in (col + 1).. Matrix raise MatrixError {
  self.row_del(row).col_del(col)
}

///|
/// Compute one matrix minor.
///
/// - Does: Takes the determinant of the submatrix formed by removing one row and one column.
/// - Input: A `Matrix`, a row index, and a column index.
/// - Returns: One symbolic `Expr`.
/// - Limits: Propagates index and determinant errors from `minor_submatrix()` and `det()`.
pub fn Matrix::minor(
  self : Matrix,
  row : Int,
  col : Int,
) -> Expr raise MatrixError {
  self.minor_submatrix(row, col).det()
}

///|
/// Compute one cofactor.
///
/// - Does: Multiplies the selected minor by the alternating sign `(-1)^(row + col)`.
/// - Input: A `Matrix`, a row index, and a column index.
/// - Returns: One symbolic `Expr`.
/// - Limits: Propagates index and determinant errors from `minor()`.
pub fn Matrix::cofactor(
  self : Matrix,
  row : Int,
  col : Int,
) -> Expr raise MatrixError {
  let sign = if (row + col) % 2 == 0 {
    @symcore.int(1)
  } else {
    expr_neg(@symcore.int(1))
  }
  expr_mul(sign, self.minor(row, col))
}

///|
/// Compute the cofactor matrix.
///
/// - Does: Replaces each entry with its cofactor.
/// - Input: A square dense `Matrix`.
/// - Returns: A dense `Matrix` of the same shape.
/// - Limits: Raises `MatrixError::NonSquareMatrixError` when the matrix is not square.
pub fn Matrix::cofactor_matrix(self : Matrix) -> Matrix raise MatrixError {
  if !self.is_square() {
    raise MatrixError::NonSquareMatrixError(
      "cofactor matrix requires a square matrix",
    )
  }
  let out : Array[Array[Expr]] = []
  for i in 0.. Matrix raise MatrixError {
  self.cofactor_matrix().transpose()
}

///|
/// Compute the inverse of a square matrix.
///
/// - Does: Chooses an exact or numeric elimination path and returns the matrix inverse.
/// - Input: A square dense `Matrix`.
/// - Returns: A dense `Matrix` whose product with the input is the identity when the input is invertible.
/// - Limits: Raises `MatrixError::NonSquareMatrixError` for non-square inputs and `MatrixError::SingularMatrixError` for singular matrices.
///
/// ```mbt check
/// test "inv computes a dense matrix inverse" {
///   let m = @symmatrices.matrix([
///     [@symcore.int(1), @symcore.int(2)],
///     [@symcore.int(3), @symcore.int(4)],
///   ])
///   inspect(m.inv().to_string(), content="Matrix([[-2, 1], [3/2, -1/2]])")
/// }
/// ```
pub fn Matrix::inv(self : Matrix) -> Matrix raise MatrixError {
  if !self.is_square() {
    raise MatrixError::NonSquareMatrixError("inverse requires a square matrix")
  }
  let n = self.rows
  if !matrix_has_floating_leaf(self.data) {
    match matrix_try_exact_inverse(self) {
      Some(inverse) => return inverse
      None => ()
    }
  }
  let numeric_prec = expr_numeric_precision(self.data)
  if matrix_has_floating_leaf(self.data) &&
    matrix_all_numeric(self.data, numeric_prec) {
    let left = clone_expr_rows(self.evalf(prec=numeric_prec).data)
    let right = clone_expr_rows(eye(n).evalf(prec=numeric_prec).data)
    for col in 0.. 0 => {
            pivot_row = candidate
            pivot_measure = Some(measure)
          }
          (Some(measure), None) => {
            pivot_row = candidate
            pivot_measure = Some(measure)
          }
          _ => ()
        }
      }
      if expr_numeric_zero(left[pivot_row][col], numeric_prec) {
        raise MatrixError::SingularMatrixError("matrix is singular")
      }
      if pivot_row != col {
        swap_dense_rows(left, pivot_row, col)
        swap_dense_rows(right, pivot_row, col)
      }
      let pivot = left[col][col]
      for j in 0.. Matrix raise MatrixError {
  if self.rows != rhs.rows {
    raise MatrixError::ShapeError(
      "solve expects rhs with the same number of rows",
    )
  }
  if !self.is_square() {
    raise MatrixError::NonSquareMatrixError("solve requires a square lhs")
  }
  if !matrix_has_floating_leaf(self.data) && !matrix_has_floating_leaf(rhs.data) {
    match matrix_try_exact_solve(self, rhs) {
      Some(solution) => return solution
      None => ()
    }
  }
  let (l, u, swaps) = self.lu()
  let permuted_rhs = clone_expr_rows(rhs.data)
  for swap in swaps {
    swap_dense_rows(permuted_rhs, swap.0, swap.1)
  }
  let y = l.lower_triangular_solve({
    rows: rhs.rows,
    cols: rhs.cols,
    data: permuted_rhs,
  })
  u.upper_triangular_solve(y)
}

///|
pub fn Matrix::rref_rhs(
  self : Matrix,
  rhs : Matrix,
) -> (Matrix, Matrix) raise MatrixError {
  if self.rows != rhs.rows {
    raise MatrixError::ShapeError(
      "rref_rhs expects rhs with the same number of rows",
    )
  }
  let augmented = hstack([self, eye(self.rows), rhs])
  let (reduced0, _) = augmented.rref()
  let reduced = reduced0.applyfunc(expr_simplify)
  let left_cols : Array[Int] = []
  let right_cols : Array[Int] = []
  for j in 0.. Array[Array[Expr]] {
  let (rref, _) = self.rref()
  let out : Array[Array[Expr]] = []
  for row in rref.data {
    if row.any(item => !expr_is_zero(item)) {
      let copied : Array[Expr] = []
      for item in row {
        copied.push(item)
      }
      out.push(copied)
    }
  }
  out
}

///|
/// Return a basis for the column space.
///
/// - Does: Selects the pivot columns from the original matrix.
/// - Input: Any dense `Matrix`.
/// - Returns: `Array[Array[Expr]]`, one basis vector per column.
/// - Limits: Uses pivot detection from `rref()`, so symbolic pivot choices follow the same heuristics.
pub fn Matrix::columnspace(self : Matrix) -> Array[Array[Expr]] {
  let (_, pivots) = self.rref()
  let out : Array[Array[Expr]] = []
  for pivot in pivots {
    let col : Array[Expr] = []
    for i in 0.. Array[Array[Expr]] {
  let (rref, pivots) = self.rref()
  let pivot_set : Map[Int, Bool] = {}
  for pivot in pivots {
    pivot_set[pivot] = true
  }
  let free_cols : Array[Int] = []
  for j in 0.. (Matrix, Matrix, Array[(Int, Int)]) raise MatrixError {
  if !self.is_square() {
    raise MatrixError::NonSquareMatrixError(
      "LU decomposition requires a square matrix",
    )
  }
  let n = self.rows
  let u = clone_expr_rows(self.data)
  let l = clone_expr_rows(eye(n).data)
  let swaps : Array[(Int, Int)] = []
  for k in 0.. Matrix raise MatrixError {
  if !self.is_square() || self.rows != rhs.rows {
    raise MatrixError::ShapeError(
      "lower_triangular_solve expects a square lhs and matching rhs rows",
    )
  }
  let numeric_prec = {
    let lhs_prec = expr_numeric_precision(self.data)
    let rhs_prec = expr_numeric_precision(rhs.data)
    if lhs_prec > rhs_prec {
      lhs_prec
    } else {
      rhs_prec
    }
  }
  let numeric_mode = (
      matrix_has_floating_leaf(self.data) || matrix_has_floating_leaf(rhs.data)
    ) &&
    matrix_all_numeric(self.data, numeric_prec) &&
    matrix_all_numeric(rhs.data, numeric_prec)
  let out = zeros(self.rows, cols=Some(rhs.cols))
  let data = clone_expr_rows(out.data)
  for j in 0.. Matrix raise MatrixError {
  if !self.is_square() || self.rows != rhs.rows {
    raise MatrixError::ShapeError(
      "upper_triangular_solve expects a square lhs and matching rhs rows",
    )
  }
  let numeric_prec = {
    let lhs_prec = expr_numeric_precision(self.data)
    let rhs_prec = expr_numeric_precision(rhs.data)
    if lhs_prec > rhs_prec {
      lhs_prec
    } else {
      rhs_prec
    }
  }
  let numeric_mode = (
      matrix_has_floating_leaf(self.data) || matrix_has_floating_leaf(rhs.data)
    ) &&
    matrix_all_numeric(self.data, numeric_prec) &&
    matrix_all_numeric(rhs.data, numeric_prec)
  let out = zeros(self.rows, cols=Some(rhs.cols))
  let data = clone_expr_rows(out.data)
  for j in 0.. Matrix raise MatrixError {
  if !self.is_square() {
    raise MatrixError::NonSquareMatrixError(
      "cholesky decomposition requires a square matrix",
    )
  }
  let n = self.rows
  let out = zeros(n, cols=Some(n))
  let l = clone_expr_rows(out.data)
  let numeric_prec = expr_numeric_precision(self.data)
  let numeric_mode = matrix_has_floating_leaf(self.data) &&
    matrix_all_numeric(self.data, numeric_prec)
  for i in 0.. (Matrix, Matrix) raise MatrixError {
  if !self.is_square() {
    raise MatrixError::NonSquareMatrixError(
      "LDL decomposition requires a square matrix",
    )
  }
  let n = self.rows
  let l0 = zeros(n, cols=Some(n))
  let d0 = zeros(n, cols=Some(n))
  let l = clone_expr_rows(l0.data)
  let d = clone_expr_rows(d0.data)
  for i in 0.. (Matrix, Matrix) raise MatrixError {
  if self.rows == 0 || self.cols == 0 {
    return (
      zeros(self.rows, cols=Some(self.cols)),
      zeros(self.cols, cols=Some(self.cols)),
    )
  }
  let q = clone_expr_rows(self.data)
  let r = Array::makei(self.cols, fn(_) {
    Array::make(self.cols, @symcore.int(0))
  })
  let ranked : Array[Int] = []
  for j in 0..