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

module Math.Projects.ChevalleyGroup.Exceptional where

import Prelude hiding ( (*>) )

import Data.List as L

-- import Math.Algebra.Field.Base
-- import Math.Algebra.Field.Extension hiding ( (<+>), (<*>) )
import Math.Core.Field
import Math.Algebra.LinearAlgebra

import Math.Algebra.Group.PermutationGroup hiding (fromList)
-- import Math.Algebra.Group.SchreierSims as SS
import Math.Algebra.Group.RandomSchreierSims as RSS

import Math.Combinatorics.FiniteGeometry (ptsAG)


-- Follows Conway's notation

-- The octonion xinf + x0 i0 + x1 i1 + ... + x6 i6
-- is represented as O [(-1,xinf),(0,x0),(1,x1),...,(6,x6)]
-- where a list element may be omitted if the coefficient is zero


newtype Octonion k = O [(Int,k)] deriving (Eq, Ord)

i0, i1, i2, i3, i4, i5, i6 :: Octonion Q
i0 = O [(0,1)]
i1 = O [(1,1)]
i2 = O [(2,1)]
i3 = O [(3,1)]
i4 = O [(4,1)]
i5 = O [(5,1)]
i6 = O [(6,1)]

fromList as = O $ filter ((/=0) . snd) $ zip [-1..6] as

toList (O xs) = toList' xs [-1..6] where
    toList' ((i,a):xs) (j:js) =
        if i == j then a : toList' xs js else 0 : toList' ((i,a):xs) js
    toList' [] (j:js) = 0 : toList' [] js
    toList' _ [] = []

expose (O ts) = ts

instance Show k => Show (Octonion k) where
    show (O []) = "0"
    show (O ts) = let c:cs = concatMap showTerm ts
                  in if c == '+' then cs else c:cs
        where showTerm (i,a) = showCoeff a ++ showImag a i
              showCoeff a = case show a of
                            "1" -> "+"
                            "-1" -> "-"
                            '-':cs -> '-':cs
                            cs -> '+':cs
              showImag a i | i == -1 = case show a of
                                       "1" -> "1"
                                       "-1" -> "1"
                                       otherwise -> ""
                           | otherwise = "i" ++ show i

instance (Ord k, Num k) => Num (Octonion k) where -- Ord k not strictly required, but keeps nf simpler
    O ts + O us = O $ nf $ ts ++ us
    negate (O ts) = O $ map (\(i,a) -> (i,-a)) ts
    O ts * O us = O $ nf [m t u | t <- ts, u <- us]
    fromInteger 0 = O []
    fromInteger n = O [(-1, fromInteger n)]

nf ts = nf' $ L.sort ts where
    nf' ((i1,a1):(i2,a2):ts) =
        if i1 == i2
        then if a1+a2 == 0 then nf' ts else nf' ((i1,a1+a2):ts)
        else (i1,a1) : nf' ((i2,a2):ts)
    nf' ts = ts

m (-1,a) (i,b) = (i,a*b)
m (i,a) (-1,b) = (i,a*b)
m (i,a) (j,b) =
    case (j-i) `mod` 7 of
    0 -> (-1,-a*b)
    1 -> ( (i+3) `mod` 7, a*b)  -- i_n+1 * i_n+2 == i_n+4
    2 -> ( (i+6) `mod` 7, a*b)  -- i_n+2 * i_n+4 == i_n+1
    3 -> ( (i+1) `mod` 7, -a*b) -- i_n+1 * i_n+4 == -i_n+2
    4 -> ( (i+5) `mod` 7, a*b)  -- i_n+4 * i_n+1 == i_n+2
    5 -> ( (i+4) `mod` 7, -a*b) -- i_n+4 * i_n+2 == -i_n+1
    6 -> ( (i+2) `mod` 7, -a*b) -- i_n+2 * i_n+1 == -i_n+4


conj (O ts) = O $ map (\(i,a) -> if i == -1 then (i,a) else (i,-a)) ts

sqnorm (O ts) = sum [a^2 | (i,a) <- ts]

instance (Ord k, Num k, Fractional k) => Fractional (Octonion k) where
    recip x = let O x' = conj x
                  xx' = sqnorm x
              in O $ map (\(i,a) -> (i,a/xx')) x'


isOrthogonal (O ts) (O us) = dot ts us == 0 where
    dot ((i,a):ts) ((j,b):us) =
        case compare i j of
        EQ -> a*b + dot ts us
        LT -> dot ts ((j,b):us)
        GT -> dot ((i,a):ts) us
    dot _ _ = 0

antiCommutes x y = x*y + y*x == 0

-- anti-commuting and being orthogonal appear to be equivalent for unit imaginary octonions,
-- provided we're not in characteristic 2


-- OCTONIONS OVER FINITE FIELDS

{-
octonions fq = map O $ octonions' [-1..6] where
    octonions' (i:is) = [if a == 0 then ts else (i,a):ts | a <- fq, ts <- octonions' is]
    octonions' [] = [[]]
-}
octonions fq = map fromList $ ptsAG 8 fq

isUnit x = sqnorm x == 1

unitImagOctonions fq = filter isUnit $ map (fromList . (0:)) $ ptsAG 7 fq


-- given the images of i0, i1, i2, return the automorphism
-- the inputs must be pure imaginary unit octonions
-- and we must have isOrthogonal i0 i1, isOrthogonal i0 i2, isOrthogonal i1 i2, and isOrthogonal (i0*i1) i2
autFrom i0' i1' i2' =
    let 0:r0 = toList i0'
        0:r1 = toList i1'
        0:r2 = toList i2'
        0:r3 = toList $ i0'*i1'
        0:r4 = toList $ i1'*i2'
        0:r5 = toList $ i0'*(i1'*i2')
        0:r6 = toList $ i0'*i2'
    in [r0,r1,r2,r3,r4,r5,r6]


x %^ g =
    let a:as = toList x
    in fromList $ a : (as <*>> g)


-- G2(3)

alpha3 = autFrom (O [(1,1::F3)]) (O [(2,1)]) (O [(3,1)])
beta3 = autFrom (O [(0,1::F3)]) (O [(2,1)]) (O [(4,1)])

gamma3s = [x | x <- unitImagOctonions f3, isOrthogonal (O [(0,1)]) x, isOrthogonal (O [(1,1)]) x, isOrthogonal (O [(3,1)]) x]

gamma3 = autFrom (O [(0,1::F3)]) (O [(1,1)]) (O [(2,1),(4,1),(5,1),(6,1)])

alpha3' = fromPairs [(x, x %^ alpha3) | x <- unitImagOctonions f3]
beta3' = fromPairs [(x, x %^ beta3) | x <- unitImagOctonions f3]
gamma3' = fromPairs [(x, x %^ gamma3) | x <- unitImagOctonions f3]

-- |Generators for G2(3), a finite simple group of order 4245696,
-- as a permutation group on the 702 unit imaginary octonions over F3
g2_3 :: [Permutation (Octonion F3)]
g2_3 = [alpha3', beta3', gamma3']
-- These three together generate a group of order 4245696, which is therefore the whole of G2(3)

-- Unit imaginary octonions form one orbit under the action of G2

-- [alpha', beta', gamma'] generate G2(3) as a permutation group on 702 points (the number of unit imaginary octonions over F3)
-- Interestingly, http://brauer.maths.qmul.ac.uk/Atlas/v3/exc/G23/ doesn't seem to have this permutation representation


-- G2(4)

alpha4 = autFrom (O [(1,1::F4)]) (O [(2,1)]) (O [(3,1)])
beta4 = autFrom (O [(0,1::F4)]) (O [(2,1)]) (O [(4,1)])

gamma4s = [x | x <- unitImagOctonions f4, isOrthogonal (O [(0,1)]) x, isOrthogonal (O [(1,1)]) x, isOrthogonal (O [(3,1)]) x]

-- gamma4 = autFrom (O [(0,1::F4)]) (O [(1,1)]) (O [(5,embed x),(6,embed $ 1+x)])
gamma4 = autFrom (O [(0,1::F4)]) (O [(1,1)]) (O [(5,a4),(6,1+a4)])

alpha4' = fromPairs [(x, x %^ alpha4) | x <- unitImagOctonions f4]
beta4' = fromPairs [(x, x %^ beta4) | x <- unitImagOctonions f4]
gamma4' = fromPairs [(x, x %^ gamma4) | x <- unitImagOctonions f4]

-- Haven't checked whether these generate whole group - can be expected to run a long time