HaskellForMaths-0.2.0: Math/Projects/MiniquaternionGeometry.hs
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
findQuadrangles' i (ls,rs)
| 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]
-}