{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

module Internal.Determinant where

import qualified Data.List as List
import Data.List.Split
import Data.Proxy
import qualified GHC.Natural as Natural
import GHC.TypeNats
import Internal.Matrix
import QLinear.Index

-- | Determinant of matrix
--
-- >>> det [matrix| 1 0; 0 1|]
-- 1
-- >>> det [matrix| 1 3; 4 2|]
-- -10
det :: Num a => Matrix n n a -> a
det :: Matrix n n a -> a
det (Matrix (1, _) [[a :: a
a]]) = a
a
det (Matrix (2, _) [[a :: a
a, b :: a
b], [c :: a
c, d :: a
d]]) = a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
c
det (Matrix (3, _) [[a :: a
a, b :: a
b, c :: a
c], [d :: a
d, e :: a
e, f :: a
f], [g :: a
g, h :: a
h, i :: a
i]]) =
  a
a a -> a -> a
forall a. Num a => a -> a -> a
* (a
e a -> a -> a
forall a. Num a => a -> a -> a
* a
i a -> a -> a
forall a. Num a => a -> a -> a
- a
f a -> a -> a
forall a. Num a => a -> a -> a
* a
h) a -> a -> a
forall a. Num a => a -> a -> a
- a
d a -> a -> a
forall a. Num a => a -> a -> a
* (a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
i a -> a -> a
forall a. Num a => a -> a -> a
- a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
h) a -> a -> a
forall a. Num a => a -> a -> a
+ a
g a -> a -> a
forall a. Num a => a -> a -> a
* (a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
f a -> a -> a
forall a. Num a => a -> a -> a
- a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
e)
det (Matrix (n :: Int
n, _) matrix :: [[a]]
matrix) = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ ((Int, Int) -> a) -> [(Int, Int)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> a
calcElem [(Int, Int)]
indices
  where
    calcElem :: (Int, Int) -> a
calcElem index :: (Int, Int)
index@(i :: Int
i, j :: Int
j) = (([a] -> Int -> a
forall a. [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1)) ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ [[a]]
matrix [[a]] -> Int -> [a]
forall a. [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1)) a -> a -> a
forall a. Num a => a -> a -> a
* [[a]] -> Int -> (Int, Int) -> a
forall a. Num a => [[a]] -> Int -> (Int, Int) -> a
unsafeAlgebraicComplement [[a]]
matrix Int
n (Int, Int)
index
    indices :: [(Int, Int)]
indices = [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Int -> [Int]
forall a. a -> [a]
repeat 1) [1 .. Int
n]

-- | Typesafe algebraic complement
--
-- To use it you have to know i and j at compile time
--
-- >>> algebraicComplement [matrix| 1 2; 3 4 |] (Index @1 @1)
-- 4
-- >>> algebraicComplement [matrix| 1 2 3; 4 5 6; 7 8 9 |] (Index @1 @1)
-- -3
algebraicComplement ::
  forall n a i j.
  (KnownNat i, KnownNat j, KnownNat n, Num a, i <= n, j <= n) =>
  Matrix n n a ->
  Index i j ->
  a
algebraicComplement :: Matrix n n a -> Index i j -> a
algebraicComplement (Matrix (n :: Int
n, _) matrix :: [[a]]
matrix) _ = [[a]] -> Int -> (Int, Int) -> a
forall a. Num a => [[a]] -> Int -> (Int, Int) -> a
unsafeAlgebraicComplement [[a]]
matrix Int
n (Int
i, Int
j)
  where
    i :: Int
i = (Natural -> Int
Natural.naturalToInt (Natural -> Int) -> Natural -> Int
forall a b. (a -> b) -> a -> b
$ Proxy i -> Natural
forall (n :: Nat) (proxy :: Nat -> *).
KnownNat n =>
proxy n -> Natural
natVal (Proxy i -> Natural) -> Proxy i -> Natural
forall a b. (a -> b) -> a -> b
$ Proxy i
forall k (t :: k). Proxy t
Proxy @i)
    j :: Int
j = (Natural -> Int
Natural.naturalToInt (Natural -> Int) -> Natural -> Int
forall a b. (a -> b) -> a -> b
$ Proxy j -> Natural
forall (n :: Nat) (proxy :: Nat -> *).
KnownNat n =>
proxy n -> Natural
natVal (Proxy j -> Natural) -> Proxy j -> Natural
forall a b. (a -> b) -> a -> b
$ Proxy j
forall k (t :: k). Proxy t
Proxy @j)

-- | Algebraic complement.
--
-- Use it if you don't know indices at compile time
--
-- >>> algebraicComplement' [matrix| 1 2; 3 4 |] (1, 1)
-- Just 4
--
-- >>> algebraicComplement' [matrix| 1 2; 3 4 |] (34, 43)
-- Nothing
--
-- >>> algebraicComplement' [matrix| 1 2 3; 4 5 6; 7 8 9 |] (1, 1)
-- Just (-3)
algebraicComplement' :: Num a => Matrix n n a -> (Int, Int) -> Maybe a
algebraicComplement' :: Matrix n n a -> (Int, Int) -> Maybe a
algebraicComplement' (Matrix (n :: Int
n, _) matrix :: [[a]]
matrix) ij :: (Int, Int)
ij@(i :: Int
i, j :: Int
j)
  | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
n Bool -> Bool -> Bool
&& Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
n = a -> Maybe a
forall a. a -> Maybe a
Just (a -> Maybe a) -> a -> Maybe a
forall a b. (a -> b) -> a -> b
$ [[a]] -> Int -> (Int, Int) -> a
forall a. Num a => [[a]] -> Int -> (Int, Int) -> a
unsafeAlgebraicComplement [[a]]
matrix Int
n (Int, Int)
ij
  | Bool
otherwise = Maybe a
forall a. Maybe a
Nothing

-- | Adjugate matrix
--
-- >>> adjugate [matrix| 1 2; 3 4|]
-- [4,-2]
-- [-3,1]
adjugate :: Num a => Matrix n n a -> Matrix n n a
adjugate :: Matrix n n a -> Matrix n n a
adjugate (Matrix size :: (Int, Int)
size@(n :: Int
n, _) matrix :: [[a]]
matrix) = (Int, Int) -> [[a]] -> Matrix n n a
forall (m :: Nat) (n :: Nat) a. (Int, Int) -> [[a]] -> Matrix m n a
Matrix (Int, Int)
size ([[a]] -> Matrix n n a) -> [[a]] -> Matrix n n a
forall a b. (a -> b) -> a -> b
$ Int -> [a] -> [[a]]
forall e. Int -> [e] -> [[e]]
chunksOf Int
n ([a] -> [[a]]) -> [a] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [a]
adj
  where
    adj :: [a]
adj = ((Int, Int) -> a) -> [(Int, Int)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ([[a]] -> Int -> (Int, Int) -> a
forall a. Num a => [[a]] -> Int -> (Int, Int) -> a
unsafeAlgebraicComplement [[a]]
matrix Int
n) [(Int
i, Int
j) | Int
j <- [1 .. Int
n], Int
i <- [1 .. Int
n]]

unsafeAlgebraicComplement :: forall a. Num a => [[a]] -> Int -> (Int, Int) -> a
unsafeAlgebraicComplement :: [[a]] -> Int -> (Int, Int) -> a
unsafeAlgebraicComplement matrix :: [[a]]
matrix n :: Int
n (i :: Int
i, j :: Int
j) = a
k a -> a -> a
forall a. Num a => a -> a -> a
* Matrix Any Any a -> a
forall a (n :: Nat). Num a => Matrix n n a -> a
det ((Int, Int) -> [[a]] -> Matrix Any Any a
forall (m :: Nat) (n :: Nat) a. (Int, Int) -> [[a]] -> Matrix m n a
Matrix (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1, Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1) [[a]]
minor)
  where
    minor :: [[a]]
minor = [[a]] -> (Int, Int) -> [[a]]
forall a. [[a]] -> (Int, Int) -> [[a]]
cross [[a]]
matrix (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1, Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1)
    k :: a
k = Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> a) -> Integer -> a
forall a b. (a -> b) -> a -> b
$ (-1) Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i)

cross :: [[a]] -> (Int, Int) -> [[a]]
cross :: [[a]] -> (Int, Int) -> [[a]]
cross matrix :: [[a]]
matrix (i :: Int
i, j :: Int
j) = [[a]]
withoutColumn
  where
    withoutLine :: [[a]]
withoutLine = Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
cut Int
i [[a]]
matrix
    withoutColumn :: [[a]]
withoutColumn = ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
cut Int
j) [[a]]
withoutLine

cut' :: Int -> [a] -> (a, [a])
cut' :: Int -> [a] -> (a, [a])
cut' 0 (x :: a
x : xs :: [a]
xs) = (a
x, [a]
xs)
cut' n :: Int
n (x :: a
x : xs :: [a]
xs) = (a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([a] -> [a]) -> (a, [a]) -> (a, [a])
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> [a] -> (a, [a])
forall a. Int -> [a] -> (a, [a])
cut' (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1) [a]
xs

cut :: Int -> [a] -> [a]
cut :: Int -> [a] -> [a]
cut 0 (x :: a
x : xs :: [a]
xs) = [a]
xs
cut n :: Int
n (x :: a
x : xs :: [a]
xs) = a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
cut (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- 1) [a]
xs