///|
pub struct MpfMatrix {
  rows : Int
  cols : Int
  data : Array[@mpf.RawMpf]
} derive(Debug, Eq)

///|
pub impl Show for MpfMatrix with fn output(self, logger) {
  logger.write_object(self.to_repr())
}

///|
fn p3_idx(cols : Int, r : Int, c : Int) -> Int {
  r * cols + c
}

///|
fn p3_copy_vec(xs : ArrayView[@mpf.RawMpf]) -> Array[@mpf.RawMpf] {
  let out : Array[@mpf.RawMpf] = []
  for x in xs {
    out.push(x)
  }
  out
}

///|
fn p3_check_dims(rows : Int, cols : Int, name : String) -> Unit raise MPError {
  if rows < 0 || cols < 0 {
    raise ValueError("\{name}: dimensions must be non-negative")
  }
}

///|
fn p3_matrix_data_with_fill(
  rows : Int,
  cols : Int,
  fill : @mpf.RawMpf,
) -> Array[@mpf.RawMpf] {
  let out : Array[@mpf.RawMpf] = []
  for _ in 0..<(rows * cols) {
    out.push(fill)
  }
  out
}

///|
fn p3_swap_rows(
  data : Array[@mpf.RawMpf],
  cols : Int,
  r0 : Int,
  r1 : Int,
) -> Unit {
  if r0 == r1 {
    return
  }
  for j in 0.. Unit {
  if c0 == c1 {
    return
  }
  for i in 0.. Int raise MPError {
  if r < 0 || r >= mat.rows || c < 0 || c >= mat.cols {
    raise ValueError("\{name}: index out of range")
  }
  p3_idx(mat.cols, r, c)
}

///|
pub fn MPContext::matrix(
  self : MPContext,
  rows : Int,
  cols : Int,
  fill? : @mpf.RawMpf = @mpf.fzero,
) -> MpfMatrix raise MPError {
  p3_check_dims(rows, cols, "matrix")
  let v = @mpf.mpf_pos(fill, self.precision(), self.round_mode())
  { rows, cols, data: p3_matrix_data_with_fill(rows, cols, v) }
}

///|
pub fn MPContext::matrix_from_rows(
  self : MPContext,
  rows_data : ArrayView[Array[@mpf.RawMpf]],
) -> MpfMatrix raise MPError {
  ignore(self)
  if rows_data.length() == 0 {
    return { rows: 0, cols: 0, data: [] }
  }
  let rows = rows_data.length()
  let cols = rows_data[0].length()
  let data : Array[@mpf.RawMpf] = []
  for i in 0.. MpfMatrix raise MPError {
  p3_check_dims(n, n, "eye")
  let out = self.matrix(n, n)
  for i in 0.. @mpf.RawMpf raise MPError {
  self.data[p3_check_index(self, r, c, "matrix.get")]
}

///|
pub fn MpfMatrix::set(
  self : MpfMatrix,
  r : Int,
  c : Int,
  v : @mpf.RawMpf,
) -> Unit raise MPError {
  self.data[p3_check_index(self, r, c, "matrix.set")] = v
}

///|
pub fn MPContext::matrix_transpose(
  self : MPContext,
  a : MpfMatrix,
) -> MpfMatrix raise MPError {
  let out = self.matrix(a.cols, a.rows)
  for i in 0.. MpfMatrix raise MPError {
  if a.rows != b.rows || a.cols != b.cols {
    raise ValueError("matrix_add: incompatible dimensions")
  }
  let p = self.p2_work_prec()
  let out = self.matrix(a.rows, a.cols)
  for i in 0..<(a.rows * a.cols) {
    out.data[i] = @mpf.mpf_add(a.data[i], b.data[i], p, @mpf.round_nearest)
  }
  out
}

///|
pub fn MPContext::matrix_sub(
  self : MPContext,
  a : MpfMatrix,
  b : MpfMatrix,
) -> MpfMatrix raise MPError {
  if a.rows != b.rows || a.cols != b.cols {
    raise ValueError("matrix_sub: incompatible dimensions")
  }
  let p = self.p2_work_prec()
  let out = self.matrix(a.rows, a.cols)
  for i in 0..<(a.rows * a.cols) {
    out.data[i] = @mpf.mpf_sub(a.data[i], b.data[i], p, @mpf.round_nearest)
  }
  out
}

///|
pub fn MPContext::matrix_mul(
  self : MPContext,
  a : MpfMatrix,
  b : MpfMatrix,
) -> MpfMatrix raise MPError {
  if a.cols != b.rows {
    raise ValueError("matrix_mul: incompatible dimensions")
  }
  let p = self.p2_work_prec()
  let out = self.matrix(a.rows, b.cols)
  for i in 0.. Array[@mpf.RawMpf] raise MPError {
  if x.length() != a.cols {
    raise ValueError("matrix_vec_mul: incompatible dimensions")
  }
  let p = self.p2_work_prec()
  let out : Array[@mpf.RawMpf] = []
  for i in 0.. @mpf.RawMpf {
  let p = self.p2_work_prec()
  let mut best = @mpf.fzero
  for i in 0.. Array[@mpf.RawMpf] raise MPError {
  if a.rows != a.cols {
    raise ValueError("lu_solve: matrix must be square")
  }
  let n = a.rows
  if b.length() != n {
    raise ValueError("lu_solve: incompatible rhs dimension")
  }
  if n == 0 {
    return []
  }
  let p = self.p2_work_prec()
  let data = p3_copy_vec(a.data)
  let rhs = p3_copy_vec(b)
  for k in 0.. @mpf.RawMpf raise MPError {
  if a.rows != a.cols {
    raise ValueError("det: matrix must be square")
  }
  let n = a.rows
  if n == 0 {
    return @mpf.fone
  }
  let p = self.p2_work_prec()
  let data = p3_copy_vec(a.data)
  let mut sign = 1
  for k in 0.. MpfMatrix raise MPError {
  if a.rows != a.cols {
    raise ValueError("inverse: matrix must be square")
  }
  let n = a.rows
  let inv = self.matrix(n, n)
  for j in 0.. Int raise MPError {
  let m = a.rows
  let n = a.cols
  if m == 0 || n == 0 {
    return 0
  }
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => {
      let s = if self.precision() > 0 { -(self.precision() / 2) } else { -40 }
      @mpf.from_man_exp(1N, s, 0, @mpf.round_down)
    }
  }
  let data = p3_copy_vec(a.data)
  let mut row = 0
  for col in 0..= m {
      break
    }
    let mut pivot = row
    let mut pivot_abs = @mpf.fzero
    for i in row..