The Analogical Engine

 
import MyPermutations
import MyMatrices
 
type Polynomial = [Int]
 
polySum :: Polynomial -> Polynomial -> Polynomial
polySum p q = if (length p >= length q) 
	then zipWith (+) p (q ++ [0,0..])
	else zipWith (+) q (p ++ [0,0..])
 
polyProduct :: Polynomial -> Polynomial -> Polynomial
polyProduct p q = map (polyProduct' p q) [0..((length p) + (length q) - 2)]
 
polyProduct' :: Polynomial -> Polynomial -> Int -> Int
polyProduct' p q k = sum [ ((p++[0,0..])!!i) * ((q++[0,0..])!!(k-i)) | 
                                         i <- [0..k] ]
 
.

 
type MatrixPoly = [[Polynomial]]
 
polyMatTranspose :: MatrixPoly -> MatrixPoly
polyMatTranspose ([]:xs) = []
polyMatTranspose xs = (map head xs):(polyMatTranspose (map tail xs))
 
polyMatProduct :: MatrixPoly -> MatrixPoly -> MatrixPoly
polyMatProduct m1 m2 = polyMatProduct' m1 (polyMatTranspose m2)
 
polyMatProduct' :: MatrixPoly -> MatrixPoly -> MatrixPoly
polyMatProduct' [] m = []
polyMatProduct' (r:rs) m = (polyFindRow' r m):(polyMatProduct' rs m)
 
polyFindRow' :: [Polynomial] -> MatrixPoly -> [Polynomial]
polyFindRow' r m = map (foldr (polySum) [0]) (map (zipWith (polyProduct) r) m)
 
polyDeterminant :: MatrixPoly -> Polynomial
polyDeterminant m = foldr (polySum) [0] (map (polyDetTerm m) 
                                           (allPermutedLists (length m)))
 
polyDetTerm :: MatrixPoly -> PermutedList -> Polynomial
polyDetTerm m p = map (\x -> (-1)^(sign p) * x) (polyDetTerm' m p)
 
polyDetTerm' :: MatrixPoly -> PermutedList -> Polynomial
polyDetTerm' _ [] = [1]
polyDetTerm' (m:ms) (p:ps) = polyProduct (m!!(p-1)) (polyDetTerm' ms ps)
 
polyMinor :: MatrixPoly -> [Int] -> [Int] -> MatrixPoly
polyMinor m rows cols = [ [ (m!!(i-1))!!(j-1) | j <- cols ] | i <- rows ] 
 
almostPolyMatInverse :: MatrixPoly -> MatrixPoly
almostPolyMatInverse m =  [ [ polyDeterminant $ polyMinor m (foo n i) (foo n j) 
                                        |  i <- [1..n] ] | j <- [1..n]]
	where n = length m
 
.

 
charpoly :: Matrix -> Polynomial
charpoly m = polyDeterminant $ matPolyAdd (matToPolyMat m) 
      (diagonalPolyMat (take (length m) (repeat [0,-1])))
 
matToPolyMat :: Matrix -> MatrixPoly
matToPolyMat = map (map (\x -> [x])) 
 
matPolyAdd :: MatrixPoly -> MatrixPoly -> MatrixPoly
matPolyAdd = zipWith (zipWith (polySum))
 
diagonalPolyMat :: [Polynomial] -> MatrixPoly
diagonalPolyMat xs = diagonalPolyMat' (length xs) 1 xs
 
diagonalPolyMat' :: Int -> Int -> [Polynomial] -> MatrixPoly
diagonalPolyMat' _ _ [] = []
diagonalPolyMat' n k (x:xs) = ((take (k-1) (repeat [0])) ++ [x] 
         ++ (take (n-k) (repeat [0]))):(diagonalPolyMat' n (k+1) xs)
 
.