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

module Math.Projects.MiniquaternionGeometry where

import qualified Data.List as L

import Math.Common.ListSet as LS

import Math.Algebra.Field.Base
import Math.Combinatorics.FiniteGeometry (pnf, ispnf, orderPGL)
import Math.Combinatorics.Graph (combinationsOf)
import Math.Combinatorics.GraphAuts
import Math.Algebra.Group.PermutationGroup hiding (order)
import qualified Math.Algebra.Group.SchreierSims as SS
import Math.Algebra.Group.RandomSchreierSims
import Math.Combinatorics.Design as D
import Math.Algebra.LinearAlgebra -- ( (<.>), (<+>) )

import Math.Projects.ChevalleyGroup.Classical

import Test.QuickCheck

-- Sources:
-- Miniquaternion Geometry, Room & Kirkpatrick
-- Survey of Non-Desarguesian Planes, Charles Weibel

-- F9, defined by adding sqrt of -1 to F3. (The Conway poly for F9 is not so convenient for us here)
data F9 = F9 F3 F3 deriving (Eq,Ord)

instance Show F9 where
show (F9 0 0) = "0"
show (F9 0 1) = "e"
show (F9 0 2) = "-e"
show (F9 1 0) = "1"
show (F9 1 1) = "1+e"
show (F9 1 2) = "1-e"
show (F9 2 0) = "-1"
show (F9 2 1) = "-1+e"
show (F9 2 2) = "-1-e"

e = F9 0 1 -- sqrt of -1

instance Num F9 where
F9 a1 b1 + F9 a2 b2 = F9 (a1+a2) (b1+b2)
F9 a1 b1 * F9 a2 b2 = F9 (a1*a2-b1*b2) (a1*b2+a2*b1)
negate (F9 a b) = F9 (negate a) (negate b)
fromInteger n = F9 (fromInteger n) 0

f9 = [F9 a b | a <- f3, b <- f3]

w = 1-e -- a primitive element - generates the multiplicative group

conj (F9 a b) = F9 a (-b)
-- This is just the Frobenius aut x -> x^3

norm (F9 a b) = a^2 + b^2
-- == x * conj x

instance Fractional F9 where
recip x@(F9 a b) = F9 (a/n) (-b/n) where n = norm x

instance FiniteField F9 where
basisFq _ = [1,e]

-- J9, or Q, defined by modifying the multiplication in F9
data J9 = J9 F9 deriving (Eq,Ord)

instance Show J9 where
show (J9 (F9 0 0)) = "0"
show (J9 (F9 0 1)) = "-j"
show (J9 (F9 0 2)) = "j"
show (J9 (F9 1 0)) = "1"
show (J9 (F9 1 1)) = "-k"
show (J9 (F9 1 2)) = "i"
show (J9 (F9 2 0)) = "-1"
show (J9 (F9 2 1)) = "-i"
show (J9 (F9 2 2)) = "k"

squaresF9 = [1,w^2,w^4,w^6] -- and 0, but not needed here

instance Num J9 where
J9 x + J9 y = J9 (x+y)
0 * _ = 0
_ * 0 = 0
J9 x * J9 y =
if y `elem` squaresF9
then J9 (x*y)
else J9 (conj x * y)
negate (J9 x) = J9 (negate x)
fromInteger n = J9 (fromInteger n)

i = J9 w
j = J9 (w^6) -- == i-1
k = J9 (w^7) -- == i+1

j9 = [J9 x | x <- f9]

-- the aut of J9 that sends i to x
autJ9 x = fromPairs [ (a+b*i, a+b*x) | a <- [0,1,-1], b <- [1,-1] ]

autA = autJ9 (-i) -- sends i -> -i
autB = autJ9 (-k) -- sends j -> -j
autC = autJ9 (-j) -- sends k -> -k

autsJ9 = [autA, autC]
-- these two auts generate the group, which is isomorphic to S3
-- indeed, the auts permute the pairs {i,-i}, {j,-j}, {k,-k}

conj' (J9 x) = J9 (conj x)
-- Note that conj' x == x .^ autB

isAut k sigma = and [sigma x + sigma y == sigma (x+y) | x <- k, y <- k]
&& and [sigma x * sigma y == sigma (x*y) | x <- k, y <- k]

isReal x = x `elem` [0,1,-1]
isComplex = not . isReal

instance Fractional J9 where
recip 0 = error "J9.recip: 0"
recip x | isReal x  = x
| otherwise = -x

instance FiniteField J9 where
basisFq _ = [1,i,j,k]
eltsFq _ = j9

-- Near fields

prop_NearField (a,b,c) =
a+(b+c) == (a+b)+c   &&  -- addition is associative
a+b == b+a           &&  -- addition is commutative
a+0 == a             &&  -- additive identity
a+(-a) == 0          &&  -- additive inverse
a*(b*c) == (a*b)*c   &&  -- multiplication is associative
a*1 == a && 1*a == a &&  -- multiplicative identity
(a+b)*c == a*c + b*c &&  -- right-distributivity
a*0 == 0

instance Arbitrary F9 where
arbitrary = do x <- arbitrary :: Gen Int
return (f9 !! (x `mod` 9))
coarbitrary = undefined -- !! only required if we want to test functions over the type

instance Arbitrary J9 where
arbitrary = do x <- arbitrary :: Gen Int
return (j9 !! (x `mod` 9))
coarbitrary = undefined -- !! only required if we want to test functions over the type

prop_NearFieldF9 (a,b,c) = prop_NearField (a,b,c) where
types = (a,b,c) :: (F9,F9,F9)

prop_NearFieldJ9 (a,b,c) = prop_NearField (a,b,c) where
types = (a,b,c) :: (J9,J9,J9)

-- PROJECTIVE PLANES

ptsPG2 r =  [ [0,0,1] ] ++ [ [0,1,x] | x <- r ] ++ [ [1,x,y] | x <- r, y <- r ]
-- if r is sorted, then so is the result

orthogonalLinesPG2 xs = L.sort [ [x | x <- xs, u <.> x == 0] | u <- xs ]

rightLinesPG2 r =
[ [0,0,1] : [ [0,1,x] | x <- r] ] ++ -- line at infinity
[ [0,0,1] : [ [1,x,y] | y <- r] | x <- r ] ++ -- vertical lines
[ [0,1,a] : [ [1,x,y] | x <- r, y <- [x*a+b] ] | a <- r, b <- r ] -- slope multiplies on the right
-- if r is sorted, then so is the result, and each line in the result

leftLinesPG2 r =
[ [0,0,1] : [ [0,1,x] | x <- r] ] ++ -- line at infinity
[ [0,0,1] : [ [1,x,y] | y <- r] | x <- r ] ++ -- vertical lines
[ [0,1,a] : [ [1,x,y] | x <- r, y <- [a*x+b] ] | a <- r, b <- r ] -- slope multiplies on the left

-- Projective plane PG2(F9)
phi = design (xs,bs) where
xs = ptsPG2 f9
bs = orthogonalLinesPG2 xs -- L.sort [ [x | x <- xs, u <.> x == 0] | u <- xs ]

-- Then the collineations of phi consist of projective transformations,
-- together with a conjugacy collineation induced by the Frobenius aut

-- alternative construction of PG2(F9) - gives same result
phi' = design (xs,bs) where
xs = ptsPG2 f9
bs = rightLinesPG2 f9

collineationsPhi = l 3 f9 ++ [fieldAut] where
D xs bs = phi
fieldAut = fromPairs [ (x , map conj x) | x <- xs ]
-- in general, this would be PSigmaL(n,Fq), whereas we want PGammaL(n,Fq). However, for F9 they coincide.
-- order 84913920

liftToGraph (D xs bs) g = fromPairs \$ [(Left x, Left (x .^ g)) | x <- xs] ++ [(Right b, Right (b -^ g)) | b <- bs]

-- This construction appears to produce a projective plane
-- (However, Room & Kirkpatrick point out that it's not really well-defined
-- - if we had chosen different quasi-homogeneous coords, we would have got different results)
-- However, it's not the same as either omega or omegaD below
omega0 = design (xs,bs) where
xs = ptsPG2 j9
bs = orthogonalLinesPG2 xs -- L.sort [ [x | x <- xs, u <.> x == 0] | u <- xs ]

-- Room & Kirkpatrick, p103
omega = design (xs,bs) where
xs = ptsPG2 j9
bs = rightLinesPG2 j9

-- another construction that produces same result (but slower)
omega2 = design (xs,bs) where
xs = ptsPG2 j9
bs =  [ l | [p,q] <- combinationsOf 2 xs, l <- [line p q], [p,q] == take 2 l]
line p q = toListSet \$ filter ispnf [(a *> p) <+> (b *> q) | a <- j9, b <- j9]

-- Room & Kirkpatrick, p107, p114
collineationsOmega =
[r]
++ [s rho sigma | rho <- j9 \\ [0], sigma <- j9 \\ [0], rho == 1 || sigma == 1]
++ [t delta epsilon | delta <- j9, epsilon <- j9, delta * epsilon == 0] -- for generators sufficient to have only one non-zero
++ [u]
++ [a lambda | lambda <- autsJ9] where
D xs bs = omega
fromMatrix m = fromPairs [ (x, pnf (x <*>> m)) | x <- xs]
r = fromMatrix [[1,0,0],[0,0,1],[0,1,0]] -- reflect in the line x = y in the affine subplane
s rho sigma = fromPairs \$ [([1,x,y], [1,x*rho,y*sigma]) | x <- j9, y <- j9]
++ [([0,1,mu],[0,1,(recip rho)*mu*sigma]) | mu <- j9]
++ [([0,0,1],[0,0,1])] -- leaves "Y" fixed
-- fromMatrix [[1,0,0],[0,rho,0],[0,0,sigma]] -- scale x,y -> rho x, sigma y
t delta epsilon = fromMatrix [[1,delta,epsilon],[0,1,0],[0,0,1]] -- translation x,y -> x+delta, y+epsilon
u = fromPairs \$ [([1,x,y], [1,x+y,x-y]) | x <- j9, y <- j9]
++ [([0,1,mu],[0,1,-mu]) | mu <- filter isComplex j9]
++ [([0,1,0],[0,1,1]), ([0,1,1],[0,1,0]), ([0,1,-1],[0,0,1]), ([0,0,1],[0,1,-1])]
-- fromMatrix [[1,0,0],[0,1,-1],[0,1,1]]
a lambda = fromPairs [ (x, map (.^ lambda) x) | x <- xs]
-- order 311040
-- (which means this is also the plane constructed in Weibel?)

-- dual plane of omega
omegaD = design (xs,bs) where
xs = ptsPG2 j9
bs = leftLinesPG2 j9

omegaD1 = D.to1n \$ dual omega
-- need proof omega /~= omegaD

omegaD2 = design (xs,bs) where
xs = ptsPG2 j9
bs =  [ l | [p,q] <- combinationsOf 2 xs, l <- [line p q], [p,q] == take 2 l]
line p q = toListSet \$ filter ispnf [(p <* a) <+> (q <* b) | a <- j9, b <- j9]

us <* x = map (*x) us

-- Room and Kirkpatrick p130
psi = design (xs,bs) where
xs = ptsPG2 j9
isReal x = all (`elem` [0,1,-1]) x
xrs = ptsPG2 [0,1,-1] -- the thirteen real points, a copy of PG2(F3) within psi
bs = toListSet [line p q | p <- xrs, q <- xs, q /= p]
line p q = L.sort \$ p : [pnf ( (p <* a) <+> q) | a <- j9]

-- Room & Kirkpatrick p137
psi2 = design (xs,bs) where
xs = ptsPG2 j9
bs = L.sort \$
[ [0,0,1] : [ [0,1,x] | x <- j9] ] ++ -- line at infinity, z=0
[ [0,0,1] : [ [1,kappa,y] | y <- j9] | kappa <- j9 ] ++ -- vertical lines x = kappa
[ [0,1,m] : [ [1,x,m*x+kappa] | x <- j9 ] | m <- [0,1,-1], kappa <- j9 ] ++ -- lines with real slope
[ [0,1,kappa] : [ [1,x,kappa*(x-r)+s] | x <- j9 ] | r <- [0,1,-1], s <- [0,1,-1], kappa <- j9 \\ [0,1,-1] ]
-- lines with complex slope

-- Room & Kirkpatrick p134-6
collineationsPsi = realProjectivities -- real transvections, generating real projectivities
++ [a lambda | lambda <- autsJ9] where
D xs bs = psi
n = 3
realTransvections = [elemTransvection n (r,c) l | r <- [1..n], c <- [1..n], r /= c, l <- [1]]
realProjectivities = [fromPairs \$ [(x, pnf (x <*>> m)) | x <- xs] | m <- realTransvections]
a lambda = fromPairs [ (x, map (.^ lambda) x) | x <- xs]
-- order 33696

-- The order of a projective plane
order (D xs bs) = length (head bs) - 1

isProjectivePlane pi = designParams pi == Just (2,(q^2+q+1,q,1))
where q = order pi

collinear (D xs bs) ys = (not . null) [b | b <- bs, all (`elem` b) ys]

-- assume p1..4 are distinct
isQuadrangle plane ps@[p1,p2,p3,p4] =
all (not . collinear plane) (combinationsOf 3 ps)

concurrent (D xs bs) ls = (not . null) [x | x <- xs, all (x `elem`) ls]

isQuadrilateral plane ls@[l1,l2,l3,l4] =
all (not . concurrent plane) (combinationsOf 3 ls)

isOval pi ps = length ps == order pi+1
&& all (not . collinear pi) (combinationsOf 3 ps)

findOvals1 pi = findOvals' 0 ([], points pi) where
n = order pi
findOvals' i (ls,rs)
| i == n+1 = [reverse ls]
| otherwise = concatMap (findOvals' (i+1))
[ (r:ls, rs') | r:rs' <- L.tails rs, all (not . collinear pi) (map (r:) (combinationsOf 2 ls)) ]
-- if we have a function to quickly generate the line through two points,
-- then we just need to see whether the third point is on it, which is much faster than testing collinearity

findQuadrangles pi = findQuadrangles' 0 ([], points pi) where
| i == 4 = [reverse ls]
| otherwise = concatMap (findQuadrangles' (i+1))
[ (r:ls, rs') | r:rs' <- L.tails rs, all (not . collinear pi) (map (r:) (combinationsOf 2 ls)) ]

findOvals pi@(D xs bs) = findOvals' 0 ([],xs) bs where
n = order pi
findOvals' i (ls,rs) bs
| i == n+1 = [reverse ls]
| otherwise = concat
[let rls = reverse (r:ls)
(notchords, chords) = L.partition (\b -> length (rls `LS.intersect` b) < 2) bs
rs'' = foldl (\\) rs' chords
-- if any line is already a chord, remove remaining points on it from further consideration
in findOvals' (i+1) (r:ls, rs'') notchords
| r:rs' <- L.tails rs]

-- Todo:
-- Code that shows that phi is Desarguesian, and omega, omegaD and psi are not
{-
-- !! NOT WORKING
-- finds apparent counterexamples in phi too
findNonDesarguesian pi@(D xs bs) =
[ [p,x,y,z,x',y',z',k,l,m] | p <- xs,
x <- xs \\ [p],
y <- xs \\ [p,x],
z <- xs \\ [p,x,y],
(not . collinear pi) [x,y,z],
x' <- line p x \\ L.sort [p,x],
y' <- line p y \\ L.sort [p,y],
z' <- line p z \\ L.sort [p,z],
(not . collinear pi) [x',y',z'],
k <- line x y `intersect` line x' y', -- will only have one element
l <- line x z `intersect` line x' z',
m <- line y z `intersect` line y' z',
(not . collinear pi) [k,l,m]  ]
where line p q = head [b | b <- bs, p `elem` b, q `elem` b]
-}
```