module CyclicCovers where

import Matrix hiding (zeroVector, zeroMatrix, idMatrix)
import RowOperations
import List

------------------------------------------------------------

idMatrix :: Int -> Matrix
idMatrix n = [ idRow i | i <- [1..n] ]
  where
  idRow i = [ delta i j | j <- [1..n] ]
  delta i j = if i == j then 1 else 0

subtractIdentityMatrix :: Matrix -> Matrix
subtractIdentityMatrix m
  = matrixSum (matrixScalarProduct (-1) (idMatrix (numRows m))) m

lambdaMatrix :: Int -> Matrix
lambdaMatrix n = [ lambdaRow i | i <- [1..n] ]
  where
  lambdaRow i = map (lambdaEntries i) [1..n]
  lambdaEntries i j
    | i == j           = 1
    | i - 1 == j       = (-1)
    | otherwise        = 0

rowTensor :: Vector -> Matrix -> Matrix
rowTensor vec m = foldr1 (zipWith (++)) blocks
  where
  blocks = zipWith (matrixScalarProduct) vec (replicate (length vec) m)

matrixTensor :: Matrix -> Matrix -> Matrix
matrixTensor m1 m2 = concat (map (\row -> rowTensor row m2) m1)

intersectionMatrix :: Int -> Int -> Matrix
intersectionMatrix p q = matrixTensor (lambdaMatrix (p-1)) (lambdaMatrix (q-1))

seifertMatrix :: Int -> Int -> Matrix
seifertMatrix p q = matrixProduct (transpose intM) (inverse intM)
  where
  intM = intersectionMatrix p q

displayCyclicGroups :: Int -> Int -> String
displayCyclicGroups order numFactors
  | order == 1 = ""
  | numFactors == 1 && order == 0 = "Z"
  | numFactors == 1                   = "Z/" ++ (show (order))
  | order == 0                    = "Z^" ++ (show numFactors)
  | otherwise = "(Z/" ++ (show (order)) ++ ")^" ++ (show numFactors)

displayAbelianGroup :: [Int] -> String
displayAbelianGroup list
  | length cleanList == 0 = "0"
  | otherwise        = drop 3 displayString
  where
  cleanList = (filter (==0) list) ++ sort ( filter (>1) (map abs list) )
  displayString = foldr (\ns string -> " + " ++ (displayBlock ns) ++ string) 
    "" (group cleanList)
  displayBlock ns = displayCyclicGroups (head ns) (length ns)

relationMatrixRange :: Int -> Int -> Int -> Int -> [Matrix] 
-- list of (seifertMatrix)^k - id for r <= k <= s relationMatrixRange p q r s
relationMatrixRange p q r s = map subtractIdentityMatrix (
  drop r ( take (s+1) (
  iterate (matrixProduct (seifertMatrix p q)) (idMatrix( (p-1)*(q-1)) ) ) ) )

homologyCalculation :: Int -> Matrix -> String
-- for a relation matrix, calculate homology and make into a string
homologyCalculation k m =  "  k = " ++ (show k) ++ "  " ++ 
  ( (displayAbelianGroup . smithDiagonalEntries) m ) ++ "\n"

homology :: Int -> Int -> Int -> Int -> IO()
homology p q r s = putStr ( header ++
  concat [ homologyCalculation (r+k) (relationMs!!k) | k <- [0..(s-r)] ] )
  where
  header = "\nFor T( " ++ (show p) ++ ", " ++ (show q) ++ " ), " ++
    "the homology groups of the k-fold cyclic covers are:\n\n"
  relationMs = relationMatrixRange p q r s