///|
fn p3_min(a : Int, b : Int) -> Int {
  if a < b {
    a
  } else {
    b
  }
}

///|
fn MPContext::p3_lin_tol(self : MPContext, prec : Int) -> @mpf.RawMpf {
  ignore(self)
  let shift = if prec > 220 { -70 } else if prec > 140 { -56 } else { -44 }
  @mpf.from_man_exp(1N, shift, 0, @mpf.round_down)
}

///|
fn p3_vec_dot(
  xs : ArrayView[@mpf.RawMpf],
  ys : ArrayView[@mpf.RawMpf],
  prec : Int,
) -> @mpf.RawMpf {
  let mut s = @mpf.fzero
  for i in 0.. @mpf.RawMpf raise MPError {
  let s = p3_vec_dot(xs, xs, prec)
  if @mpf.mpf_sign(s) < 0 {
    raise DomainError("linear algebra: negative norm square")
  }
  @mpf.mpf_sqrt(s, prec, @mpf.round_nearest) catch {
    err => raise from_mpf_error(err)
  }
}

///|
fn p3_extract_col(a : MpfMatrix, col : Int) -> Array[@mpf.RawMpf] {
  let out : Array[@mpf.RawMpf] = []
  for i in 0.. Unit {
  for i in 0.. Array[@mpf.RawMpf] {
  let v : Array[@mpf.RawMpf] = []
  for i in 0.. Array[@mpf.RawMpf] {
  let out : Array[@mpf.RawMpf] = []
  for x in v {
    out.push(@mpf.mpf_mul(s, x, prec, @mpf.round_nearest))
  }
  out
}

///|
fn p3_vec_sub_scaled(
  v : Array[@mpf.RawMpf],
  q : ArrayView[@mpf.RawMpf],
  c : @mpf.RawMpf,
  prec : Int,
) -> Unit {
  for i in 0.. (@mpf.RawMpf, Array[@mpf.RawMpf]) raise MPError {
  let nrm = p3_vec_norm2(v, prec)
  if @mpf.mpf_le(nrm, tol) {
    let z : Array[@mpf.RawMpf] = []
    for _ in 0.. @mpf.RawMpf {
  let p = self.p2_work_prec()
  let mut best = @mpf.fzero
  for j in 0.. @mpf.RawMpf raise MPError {
  let p = self.p2_work_prec()
  let mut ss = @mpf.fzero
  for x in a.data {
    ss = @mpf.mpf_add(
      ss,
      @mpf.mpf_mul(x, x, p, @mpf.round_nearest),
      p,
      @mpf.round_nearest,
    )
  }
  let nrm = @mpf.mpf_sqrt(ss, p, @mpf.round_nearest) catch {
    err => raise from_mpf_error(err)
  }
  @mpf.mpf_pos(nrm, self.precision(), self.round_mode())
}

///|
pub fn MPContext::qr(
  self : MPContext,
  a : MpfMatrix,
  full? : Bool = false,
  tol? : @mpf.RawMpf,
) -> (MpfMatrix, MpfMatrix) raise MPError {
  let m = a.rows
  let n = a.cols
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => self.p3_lin_tol(p)
  }
  let k = p3_min(m, n)
  if k == 0 {
    return if full {
      (self.eye(m), self.matrix(m, n))
    } else {
      (self.matrix(m, 0), self.matrix(0, n))
    }
  }
  let qcols : Array[Array[@mpf.RawMpf]] = []
  let r = self.matrix(k, n)
  for j in 0.. (Array[@mpf.RawMpf], @mpf.RawMpf) raise MPError {
  let m = a.rows
  let n = a.cols
  if b.length() != m {
    raise ValueError("qr_solve: incompatible rhs dimension")
  }
  if m < n {
    raise ValueError("qr_solve: underdetermined systems not supported")
  }
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => self.p3_lin_tol(p)
  }
  let (q, r) = self.qr(a, tol=tol_abs)
  let y = p3_matrix_data_with_fill(1, n, @mpf.fzero)
  for i in 0.. raise from_mpf_error(err)
  }
  (x, @mpf.mpf_pos(rn, self.precision(), self.round_mode()))
}

///|
pub fn MPContext::cholesky(
  self : MPContext,
  a : MpfMatrix,
  tol? : @mpf.RawMpf,
) -> MpfMatrix raise MPError {
  if a.rows != a.cols {
    raise ValueError("cholesky: matrix must be square")
  }
  let n = a.rows
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => self.p3_lin_tol(p)
  }
  let l = self.matrix(n, n)
  for i in 0.. raise from_mpf_error(err)
        }
      } else {
        let ljj = l.data[p3_idx(l.cols, j, j)]
        if @mpf.mpf_le(@mpf.mpf_abs(ljj, p, @mpf.round_nearest), tol_abs) {
          raise DomainError("cholesky: matrix is singular")
        }
        l.data[p3_idx(l.cols, i, j)] = p2_mpf_div(s, ljj, p, @mpf.round_nearest)
      }
    }
  }
  for i in 0.. Array[@mpf.RawMpf] raise MPError {
  if a.rows != a.cols {
    raise ValueError("cholesky_solve: matrix must be square")
  }
  let n = a.rows
  if b.length() != n {
    raise ValueError("cholesky_solve: incompatible rhs dimension")
  }
  let p = self.p2_work_prec()
  let l = match tol {
    Some(t) => self.cholesky(a, tol=t)
    None => self.cholesky(a)
  }
  let y = p3_matrix_data_with_fill(1, n, @mpf.fzero)
  for i in 0.. (MpfMatrix, MpfMatrix) raise MPError {
  if a.rows != a.cols {
    raise ValueError("hessenberg: matrix must be square")
  }
  let n = a.rows
  if n <= 2 {
    return (self.eye(n), a)
  }
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => self.p3_lin_tol(p)
  }
  let h = p3_copy_vec(a.data)
  let q0 = self.eye(n)
  let q = p3_copy_vec(q0.data)
  for k in 0..<(n - 2) {
    let m = n - k - 1
    let x : Array[@mpf.RawMpf] = []
    for i in 0..= 0 {
      @mpf.mpf_neg(normx, p, @mpf.round_nearest)
    } else {
      normx
    }
    let u = p3_copy_vec(x)
    u[0] = @mpf.mpf_sub(u[0], alpha, p, @mpf.round_nearest)
    let (_, u_norm) = p3_vec_normalize(u, tol_abs, p)
    if @mpf.mpf_le(p3_vec_norm2(u_norm, p), tol_abs) {
      continue
    }
    for j in k.. (MpfMatrix, MpfMatrix) raise MPError {
  if a.rows != a.cols {
    raise ValueError("schur: matrix must be square")
  }
  let n = a.rows
  if n <= 1 {
    return (self.eye(n), a)
  }
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => self.p3_lin_tol(p)
  }
  let (q0, h0) = self.hessenberg(a, tol=tol_abs)
  let mut qtot = q0
  let mut t = h0
  let steps = if max_steps > 0 { max_steps } else { n * n * 96 }
  for _ in 0.. Array[@mpf.RawMpf] raise MPError {
  let n = a.rows
  let x : Array[@mpf.RawMpf] = []
  for _ in 0.. mu = @mpf.mpf_mul_int(mu, 10, prec, @mpf.round_nearest)
    } noraise {
      y => {
        let (_, next) = p3_vec_normalize(y, tol, prec)
        let mut diff = @mpf.fzero
        for i in 0.. (Array[@mpf.RawMpf], MpfMatrix) raise MPError {
  if a.rows != a.cols {
    raise ValueError("eig: matrix must be square")
  }
  let n = a.rows
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => self.p3_lin_tol(p)
  }
  let (_q, t) = self.schur(a, tol=tol_abs, max_steps~)
  for i in 1.. (Array[@mpf.RawMpf], MpfMatrix) raise MPError {
  match tol {
    Some(t) => self.eigsy(a, tol=t, max_steps~)
    None => self.eigsy(a, max_steps~)
  }
}

///|
pub fn MPContext::svd(
  self : MPContext,
  a : MpfMatrix,
  full? : Bool = false,
  tol? : @mpf.RawMpf,
) -> (MpfMatrix, Array[@mpf.RawMpf], MpfMatrix) raise MPError {
  ignore(full)
  let m = a.rows
  let n = a.cols
  let p = self.p2_work_prec()
  let tol_abs = match tol {
    Some(v) => @mpf.mpf_abs(v, p, @mpf.round_nearest)
    None => self.p3_lin_tol(p)
  }
  let k = p3_min(m, n)
  if k == 0 {
    return (self.matrix(m, 0), [], self.matrix(0, n))
  }
  if m >= n {
    let at = self.matrix_transpose(a)
    let ata = self.matrix_mul(at, a)
    let (evals, evecs) = self.eigsy(ata, tol=tol_abs)
    let u = self.matrix(m, k)
    let v = self.matrix(k, n)
    let sigmas : Array[@mpf.RawMpf] = []
    for t in 0.. raise from_mpf_error(err)
      }
      sigmas.push(@mpf.mpf_pos(sigma, self.precision(), self.round_mode()))
      let vi = p3_extract_col(evecs, idx)
      for j in 0.. raise from_mpf_error(err)
      }
      sigmas.push(@mpf.mpf_pos(sigma, self.precision(), self.round_mode()))
      let ui = p3_extract_col(uvecs, idx)
      p3_set_col(u, t, ui)
      if @mpf.mpf_le(sigma, tol_abs) {
        continue
      }
      let atui = self.matrix_vec_mul(at, ui)
      let inv = p2_mpf_div(@mpf.fone, sigma, p, @mpf.round_nearest)
      let vi = p3_vec_scale(atui, inv, p)
      for j in 0.. @mpf.RawMpf raise MPError {
  if a.rows != a.cols {
    raise ValueError("cond: matrix must be square")
  }
  let key = norm.to_lower()
  let an = match key {
    "inf" => self.matrix_norm_inf(a)
    "1" => self.matrix_norm_1(a)
    "f" => self.matrix_norm_fro(a)
    "fro" => self.matrix_norm_fro(a)
    _ => raise ValueError("cond: unknown norm '\{norm}'")
  }
  let ainv = self.inverse(a)
  let bn = match key {
    "inf" => self.matrix_norm_inf(ainv)
    "1" => self.matrix_norm_1(ainv)
    "f" => self.matrix_norm_fro(ainv)
    "fro" => self.matrix_norm_fro(ainv)
    _ => raise ValueError("cond: unknown norm '\{norm}'")
  }
  @mpf.mpf_mul(an, bn, self.p2_work_prec(), @mpf.round_nearest)
}