-- Copyright (c) David Amos, 2008. All rights reserved.

{-# OPTIONS_GHC -XFlexibleInstances #-}

module Math.Projects.KnotTheory.IwahoriHecke where

-- import qualified Data.Map as M

import Math.Algebra.Field.Base
import Math.Algebra.NonCommutative.NCPoly as NP hiding (x,y,z)
import Math.Algebra.NonCommutative.GSBasis

import Math.Projects.KnotTheory.LaurentMPoly as LP hiding (z)
import Math.Projects.KnotTheory.Braid


-- IWAHORI-HECKE ALGEBRAS

data IwahoriHeckeGens = T Int deriving (Eq,Ord)

instance Show IwahoriHeckeGens where
    show (T i) = 't': show i

t_ i = NP [(M [T i], 1)] :: NPoly LPQ IwahoriHeckeGens

t1 = t_ 1
t2 = t_ 2
t3 = t_ 3
t4 = t_ 4

q = LP.var "q" -- :: LPQ
z = LP.var "z" -- :: LPQ

q' = NP.inject q :: NPoly LPQ IwahoriHeckeGens
z' = NP.inject z :: NPoly LPQ IwahoriHeckeGens

-- inverses
instance Invertible (NPoly LPQ IwahoriHeckeGens) where
    inv (NP [(M [T i], 1)]) = (t_ i - z') / q'

-- x ^- n = inv x ^ n

-- Iwahori-Hecke algebra Hn(q,z), generated by n-1 elts t1..t_n-1, together with relations
ihRelations n =
    [t_ i * t_ j - t_ j * t_ i | i <- [1..n-1], j <- [i+2..n-1] ] ++
    [t_ i * t_ j * t_ i - t_ j * t_ i * t_ j | i <- [1..n-1], j <- [1..n-1], abs (i-j) == 1 ] ++
    [(t_ i)^2 - z' * t_ i - q' | i <- [1..n-1] ]

-- given an elt of the Temperley-Lieb algebra, return the dimension it's defined over (ie the number of points)
dimIH (NP ts) = 1 + maximum (0 : [i | (M bs,c) <- ts, T i <- bs])

-- Reduce to normal form
ihnf f = f %% (gb $ ihRelations $ dimIH f)

-- Monomial basis for Iwahori-Hecke algebra (as quotient of free algebra by Iwahori-Hecke relations)
ihBasis n = mbasisQA [t_ i | i <- [1..n-1]] (gb $ ihRelations n)


-- OCNEANU TRACE

-- Trace function on single terms
tau' 1 (1,c) = c
tau' 2 (1,c) = c * (1-q)/z
tau' 2 (m,c) = c
tau' n (m,c) = case m `divM` M [T (n-1)] of
               Just (l,r) -> tau (n-1) $ NP [(l*r,c)] -- have to call tau not tau', because l*r might be reducible
               Nothing -> tau' (n-1) (m,c*(1-q)/z)

-- Trace function on polynomials, by linearity
-- Given an elt of Iwahori-Hecke algebra, returns an elt of Q[q,z]
tau n f | dimIH f <= n = let NP ts = ihnf f in sum [tau' n t | t <- ts]


fromBraid f = ihnf (NP.subst skeinRelations f) where
    skeinRelations = concat [ [(s_ i, t_ i), (s_ (-i), (t_ i - z') / q')] | i <- [1..] ]


-- HOMFLY polynomial, also called Jones-Conway polynomial
-- Kassel, Turaev version, as poly in x,y
-- Satisfies skein relation x P_L_+(x,y) - x^-1 P_L_-(x,y) == y P_L_0(x,y)
-- n is the index of the Iwahori-Hecke algebra we're working in (ie the number of strings in the braid)
-- f is the braid expressed in the Iwahori-Hecke generators ti
homfly n f = LP.subst [(q,1/x^2),(z,y/x)] $ tau n $ fromBraid f

i = LP.var "i" :: LPQ
l = LP.var "l" :: LPQ
m = LP.var "m" :: LPQ

-- HOMFLY polynomial (for braid f over n strings)
-- Lickorish version, as poly in l,m
-- Satisfies skein relation l P_L_+(l,m) + l^-1 P_L_-(l,m) + m P_L_0(l,m) == 0
homfly' n f =
    let f' = LP.subst [(x,i^3*l),(y,i*m)] (homfly n f)
    in reduceLP f' (i^2+1)

-- Closer to Thistlethwaite notation
homfly'' n f = sum $ zipWith (*) (map LP.inject $ coeffs (m^2) (homfly' n f)) (iterate (*(m'^2)) 1)
    where m' = LP.var "m" :: LaurentMPoly LPQ
    -- where m' = LP [(LM $ M.singleton "m" 1, 1)] :: LaurentMPoly LPQ

-- express an lpoly as a upoly over one of its variables
coeffs v 0 = []
coeffs v f = let (f',c) = quotRemLP f v in c : coeffs v f'


-- Jones polynomial (for braid f over m strings)
-- from the HOMFLY polynomial via the substitution x -> 1/t, y -> t^1/2 - t^-1/2
-- The reason the code is so complicated is that y == t^1/2 - t^-1/2 can appear in the denominator,
-- and we have to cancel it out with the numerator by hand because our code can't do it for us
jones' m f = let f' = homfly m f
                 n = d*f'
                 d = denominatorLP f'
                 subs = [(x,1/t),(y,t^^^(1/2)-1/t^^^(1/2))]
                 -- subs = [(x,1/t^2),(y,t-1/t)]
                 n' = LP.subst subs n
                 d' = LP.subst subs d
                 nn = nd*n'
                 nd = denominatorLP n'
                 dn = dd*d'
                 dd = denominatorLP d'
                 (q,r) = quotRemLP nn dn
             -- in if r == 0 then halfExponents' (dd/nd * q) else error ""
             in if r == 0 then (dd/nd * q) else error ""

-- Alexander polynomial is the substitution x -> 1, y -> t^1/2 - t^-1/2