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

    {-# LANGUAGE FlexibleInstances #-}

    module Math.Projects.KnotTheory.TemperleyLieb where

    import Data.List ( (\\) )

    import Math.Algebra.Field.Base
    import Math.Algebra.NonCommutative.NCPoly as NP
    import Math.Algebra.NonCommutative.GSBasis

    import Math.Projects.KnotTheory.LaurentMPoly as LP
    import Math.Projects.KnotTheory.Braid

    -- TEMPERLEY-LIEB ALGEBRAS

    data TemperleyLiebGens = E Int deriving (Eq,Ord)

    instance Show TemperleyLiebGens where
        show (E i) = 'e': show i

    e_ i = NP [(M [E i], 1)] :: NPoly LPQ TemperleyLiebGens

    -- d is the value of a closed loop
    d = LP.var "d"
    d' = NP.inject d :: NPoly LPQ TemperleyLiebGens

    e1 = e_ 1
    e2 = e_ 2
    e3 = e_ 3
    e4 = e_ 4

    -- Temperley-Lieb algebra An(d), generated by n-1 elts e1..e_n-1, together with relations
    tlRelations n =
        [e_ i * e_ j - e_ j * e_ i | i <- [1..n-1], j <- [i+2..n-1] ] ++
        [e_ i * e_ j * e_ i - e_ i | i <- [1..n-1], j <- [1..n-1], abs (i-j) == 1 ] ++
        [(e_ i)^2 - d' * e_ i | i <- [1..n-1] ]

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

    -- Reduce to normal form
    tlnf f = f %% (gb $ tlRelations $ dimTL f)

    -- Monomial basis for Temperley-Lieb algebra (as quotient of free algebra by Temperley-Lieb relations)
    tlBasis n = mbasisQA [e_ i | i <- [1..n-1]] (gb $ tlRelations n)

    -- trace function
    -- the trace of an elt is d^k, where k is the number of loops in its closure (ie join the top and bottom of the diagram to make an annulus)
    -- this is clearly the same as the number of cycles of the elt when thought of as an elt of Sn, with ei mapped to the transposition (i i+1)
    tr' n (M g) = d ^ ( -1 + length (orbits g [1..n]) ) where
        image i [] = i
        image i (E j : es) | i == j    = image (i+1) es
                           | i == j+1  = image (i-1) es
                           | otherwise = image i es
        orbits g [] = []
        orbits g (i:is) = let i' = orbit i [] in i' : orbits g ((i:is) \\ i')
        orbit j js = let j' = image j g in if j' `elem` (j:js) then reverse (j:js) else orbit j' (j:js)
    -- Note, some authors define the trace so that tr 1 == 1.
    -- That is the same as this trace except for a factor of d^(n-1)

    tr n f@(NP ts) = sum [c * tr' n m | (m,c) <- ts]

    -- JONES POLYNOMIAL

    a = LP.var "a"
    a' = NP.inject a :: NPoly LPQ TemperleyLiebGens

    -- Convert a braid to Temperley-Lieb algebra using Skein relation
    fromBraid f = tlnf (NP.subst skeinRelations f) where
        skeinRelations = concat [ [(s_ i, 1/a' * e_ i + a'), (s_ (-i), a' * e_ i + 1/a')] | i <- [1..] ]

    -- Jones polynomial
    -- n the number of strings, f the braid
    jones n f = let kauffman = LP.subst [(d, - a^2 - 1/a^2)] $ tr n (fromBraid f)
                    j = (-a)^^(-3 * writhe f) * kauffman
                -- in halfExponents $ halfExponents $ LP.subst [(a,1/t)] j
                -- in quarterExponents' $ LP.subst [(a,1/t)] j
                in LP.subst [(a,1/t^^^(1/4))] j