// Thomas algorithm for tridiagonal systems and Householder QR decomposition.

///|
/// Thomas algorithm (Tridiagonal Matrix Algorithm - TDMA) for solving A * x = d
pub fn solve_tridiagonal_sys(
  a : Array[Double], // Lower diagonal (length n-1)
  b : Array[Double], // Main diagonal (length n)
  c : Array[Double], // Upper diagonal (length n-1)
  d : Array[Double], // Right-hand side (length n)
) -> Array[Double] {
  let n = b.length()
  if n == 0 {
    return []
  }
  let c_prime = Array::make(n, 0.0)
  let d_prime = Array::make(n, 0.0)
  let x = Array::make(n, 0.0)

  c_prime[0] = c[0] / b[0]
  d_prime[0] = d[0] / b[0]

  for i = 1; i < n; i = i + 1 {
    let m = b[i] - a[i - 1] * c_prime[i - 1]
    if i < n - 1 {
      c_prime[i] = c[i] / m
    }
    d_prime[i] = (d[i] - a[i - 1] * d_prime[i - 1]) / m
  }

  x[n - 1] = d_prime[n - 1]
  for i = n - 2; i >= 0; i = i - 1 {
    x[i] = d_prime[i] - c_prime[i] * x[i + 1]
  }

  x
}