def repeatedSquaring {a} (fn: a -> a -> a) (x: a) (n: Integer) : a :=
  match n as integer with
    | #1      -> x
    | ?isEven ->
        let y := repeatedSquaring fn x (quotient n 2)
         in fn y y
    | ?isOdd  ->
        let y := repeatedSquaring fn x (quotient n 2)
         in fn (fn y y) x


inductive pattern Integer :=
  | o
  | s Integer

def nat : Matcher Integer :=
  matcher
    | o as () with
      | 0   -> [()]
      | _   -> []
    | s $ as nat with
      | $tgt ->
          match compare tgt 0 as ordering with
            | greater -> [tgt - 1]
            | _       -> []
    | #$n as () with
      | $tgt -> if tgt = n then [()] else []
    | $ as (something) with
      | $tgt -> [tgt]


--
-- Eigenvalues and eigenvectors
--
def M.eigenvalues {Num a} (m: Matrix a) : [a] :=
  let (e1, e2) := qF (M.det (T.- m (scalarToTensor x [2, 2]))) x
   in [e1, e2]

def M.eigenvectors {Num a} (m: Matrix a) : [(a, Vector a)] :=
  let (e1, e2) := qF (M.det (T.- m (scalarToTensor x [2, 2]))) x
   in [ (e1, clearIndex (T.- m (scalarToTensor e1 [2, 2]))_i_1)
      , (e2, clearIndex (T.- m (scalarToTensor e2 [2, 2]))_i_1) ]

--
-- LU decomposition
--
def M.LU {Num a} (x: Matrix a) : (Matrix a, Matrix a) :=
  match tensorShape x as list integer with
    | [#2, #2] ->
      let L := generateTensor
                 (\[i, j] -> match compare i j as ordering with
                   | less -> 0
                   | equal -> 1
                   | greater -> b_i_j)
                 [2, 2]
          U := generateTensor
                 (\[i, j] -> match compare i j as ordering with
                   | greater -> 0
                   | _ -> c_i_j)
                 [2, 2]
          m := M.* L U
          ret := solve
                   [ (m_1_1, x_1_1, c_1_1)
                   , (m_1_2, x_1_2, c_1_2)
                   , (m_2_1, x_2_1, b_2_1)
                   , (m_2_2, x_2_2, c_2_2) ]
       in (substitute ret L, substitute ret U)
    | [#3, #3] ->
      let L := generateTensor
                 (\[i, j] -> match compare i j as ordering with
                   | less -> 0
                   | equal -> 1
                   | greater -> b_i_j)
                 [3, 3]
          U := generateTensor
                 (\[i, j] -> match compare i j as ordering with
                   | greater -> 0
                   | _ -> c_i_j)
                 [3, 3]
          m := M.* L U
          ret := solve
                   [ (m_1_1, x_1_1, c_1_1)
                   , (m_1_2, x_1_2, c_1_2)
                   , (m_1_3, x_1_3, c_1_3)
                   , (m_2_1, x_2_1, b_2_1)
                   , (m_2_2, x_2_2, c_2_2)
                   , (m_2_3, x_2_3, c_2_3)
                   , (m_3_1, x_3_1, b_3_1)
                   , (m_3_2, x_3_2, b_3_2)
                   , (m_3_3, x_3_3, c_3_3) ]
       in (substitute ret L, substitute ret U)
    | _ -> undefined

--
-- Utility
--
def generateMatrixFromQuadraticExpr {a} (f: a) (xs: [a]) : Matrix a :=
  generateTensor
    (\[i, j] -> coefficient2 f (nth i xs) (nth j xs))
    [length xs, length xs]