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

    {-# LANGUAGE FlexibleInstances #-}

    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