///|
pub(all) struct MatrixVector2 {
  x : Double
  y : Double
} derive(Eq, Debug, ToJson)

///|
pub fn MatrixVector2::new(x : Double, y : Double) -> MatrixVector2 {
  { x, y }
}

///|
pub fn MatrixVector2::add(
  a : MatrixVector2,
  b : MatrixVector2,
) -> MatrixVector2 {
  { x: a.x + b.x, y: a.y + b.y }
}

///|
pub fn MatrixVector2::scale(
  a : MatrixVector2,
  factor : Double,
) -> MatrixVector2 {
  { x: a.x * factor, y: a.y * factor }
}

///|
pub fn MatrixVector2::dot(a : MatrixVector2, b : MatrixVector2) -> Double {
  a.x * b.x + a.y * b.y
}

///|
pub fn MatrixVector2::norm(a : MatrixVector2) -> Double {
  a.dot(a).sqrt()
}

///|
pub fn MatrixVector2::distance(a : MatrixVector2, b : MatrixVector2) -> Double {
  a.add(b.scale(-1.0)).norm()
}

///|
pub(all) struct Matrix2 {
  a11 : Double
  a12 : Double
  a21 : Double
  a22 : Double
} derive(Debug, ToJson)

///|
pub fn Matrix2::new(
  a11 : Double,
  a12 : Double,
  a21 : Double,
  a22 : Double,
) -> Matrix2 {
  { a11, a12, a21, a22 }
}

///|
pub fn Matrix2::identity() -> Matrix2 {
  { a11: 1.0, a12: 0.0, a21: 0.0, a22: 1.0 }
}

///|
pub fn Matrix2::zero() -> Matrix2 {
  { a11: 0.0, a12: 0.0, a21: 0.0, a22: 0.0 }
}

///|
pub fn Matrix2::multiply_vector(
  matrix : Matrix2,
  vector : MatrixVector2,
) -> MatrixVector2 {
  {
    x: matrix.a11 * vector.x + matrix.a12 * vector.y,
    y: matrix.a21 * vector.x + matrix.a22 * vector.y,
  }
}

///|
pub fn Matrix2::add(a : Matrix2, b : Matrix2) -> Matrix2 {
  {
    a11: a.a11 + b.a11,
    a12: a.a12 + b.a12,
    a21: a.a21 + b.a21,
    a22: a.a22 + b.a22,
  }
}

///|
pub fn Matrix2::scale(a : Matrix2, factor : Double) -> Matrix2 {
  {
    a11: a.a11 * factor,
    a12: a.a12 * factor,
    a21: a.a21 * factor,
    a22: a.a22 * factor,
  }
}

///|
pub fn Matrix2::multiply(a : Matrix2, b : Matrix2) -> Matrix2 {
  {
    a11: a.a11 * b.a11 + a.a12 * b.a21,
    a12: a.a11 * b.a12 + a.a12 * b.a22,
    a21: a.a21 * b.a11 + a.a22 * b.a21,
    a22: a.a21 * b.a12 + a.a22 * b.a22,
  }
}

///|
pub fn Matrix2::determinant(matrix : Matrix2) -> Double {
  matrix.a11 * matrix.a22 - matrix.a12 * matrix.a21
}

///|
pub fn Matrix2::transpose(matrix : Matrix2) -> Matrix2 {
  { a11: matrix.a11, a12: matrix.a21, a21: matrix.a12, a22: matrix.a22 }
}

///|
pub fn Matrix2::inverse(matrix : Matrix2) -> Matrix2 {
  let determinant = matrix.determinant()
  if determinant == 0.0 {
    Matrix2::zero()
  } else {
    Matrix2::new(
      matrix.a22 / determinant,
      -matrix.a12 / determinant,
      -matrix.a21 / determinant,
      matrix.a11 / determinant,
    )
  }
}

///|
pub fn outer_product(first : MatrixVector2, second : MatrixVector2) -> Matrix2 {
  Matrix2::new(
    first.x * second.x,
    first.x * second.y,
    first.y * second.x,
    first.y * second.y,
  )
}

///|
pub fn solve_2x2(matrix : Matrix2, rhs : MatrixVector2) -> MatrixVector2 {
  matrix.inverse().multiply_vector(rhs)
}

///|
pub fn matrix_trace(matrix : Matrix2) -> Double {
  matrix.a11 + matrix.a22
}

///|
pub fn matrix_frobenius_norm(matrix : Matrix2) -> Double {
  (matrix.a11 * matrix.a11 +
  matrix.a12 * matrix.a12 +
  matrix.a21 * matrix.a21 +
  matrix.a22 * matrix.a22).sqrt()
}