-- smallGroup(162,19)

restart
-- the generators from GroupNames
a = matrix{
 {1,0,0,0,0,0},
 {0,1,0,0,0,0},
 {0,0,-1,1,0,0},
 {0,0,-1,0,0,0},
 {0,0,0,0,-1,1},
 {0,0,0,0,-1,0}
 }

b = matrix{
    {-1,1,0,0,0,0},
    {-1,0,0,0,0,0},
    {0,0,1,0,0,0},
    {0,0,0,1,0,0},
    {0,0,0,0,0,-1},
    {0,0,0,0,1,-1}
    }

c = matrix{
    {0,-1,0,0,0,0},
    {1,-1,0,0,0,0},
    {0,0,0,-1,0,0},
    {0,0,1,-1,0,0},
    {0,0,0,0,0,-1},
    {0,0,0,0,1,-1}
    }

d = matrix{
    {0,0,1,0,0,0},
    {0,0,0,1,0,0},
    {0,0,0,0,1,0},
    {0,0,0,0,0,1},
    {1,0,0,0,0,0},
    {0,1,0,0,0,0}
    }

e = matrix{
    {1,-1,0,0,0,0},
    {0,-1,0,0,0,0},
    {0,0,0,0,1,-1},
    {0,0,0,0,0,-1},
    {0,0,1,-1,0,0},
    {0,0,0,-1,0,0}
    }

-- the identity
one = id_(ZZ^6)

-- the relations from GroupNames
relationsG = (aaa,bbb,ccc,ddd,eee,one) -> (
    aaa^3 == one,
    bbb^3 == one,
    ccc^3 == one,
    ddd^3 == one,
    eee^2 == one,
    aaa*bbb == bbb*aaa,
    aaa*ccc == ccc*aaa,    
    bbb*ccc == ccc*bbb,
    ddd*ccc == ccc*ddd,
    ddd*aaa*ddd^-1==aaa*bbb,
    eee*aaa*eee==aaa^-1,
    ddd*bbb*ddd^-1==bbb*ccc,
    eee*bbb*eee==bbb*ccc^-1, 
    eee*ccc*eee==ccc^-1
    ) 

-- TEST: do the generators satisfy the relations?
assert (relationsG(a,b,c,d,e,one) == (14:true))

-- primes with a third root of unity
primes = select(100,isPrime)
select(primes,pp->3==#select(pp,i->(i^3 -1)% pp == 0))
-- {7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97}


-- a basis of eigenvectors
T = (
    syz(a-one||b-xi*one)    
    |syz(a-one||b-xi^2*one)
    |syz(a-xi*one||b-one)
    |syz(a-xi*one||b-xi^2*one)
    |syz(a-xi^2*one||b-one)
    |syz(a-xi^2*one||b-xi*one)
    )

-- apply base change to generators (diagonalizing element of C3^3)
A = T^-1*a*T
B = T^-1*b*T
C = T^-1*c*T
D = T^-1*d*T
E = T^-1*e*T
One = T^-1*one*T

-- TEST: do the generators still satisfy the relations?
assert (relationsG(A,B,C,D,E,One) == (14:true))

-- the exterior algebra
EE = K[ee_0..ee_5,SkewCommutative=>true]

-- the operation on E
opE = (w,g) -> sub(w,(vars EE)*g)

-- TEST: does the operation give back the generators?
assert (sub(diff(transpose vars EE,opE(vars EE,A)),K) == A)
assert (sub(diff(transpose vars EE,opE(vars EE,B)),K) == B)
assert (sub(diff(transpose vars EE,opE(vars EE,C)),K) == C)
assert (sub(diff(transpose vars EE,opE(vars EE,D)),K) == D)
assert (sub(diff(transpose vars EE,opE(vars EE,E)),K) == E)

-- the 2 forms
EE2 = super basis(2,EE)

-- is something G-invariant?
isInvariant = (w) -> (
    (opE(w,A) == w) and
    (opE(w,B) == w) and
    (opE(w,C) == w) and
    (opE(w,D) == w) and
    (opE(w,E) == w) 
      )

-- the weight with respect to C3^2 = <a,b,c>
weight = (w) -> ({
    coefficient(w,opE(w,A)),
    coefficient(w,opE(w,B)),
    coefficient(w,opE(w,C))
    })

-- weights
apply(flatten entries EE2,weight)

-- trivial weight vectors
w111 = matrix{{ee_0*ee_1,ee_2*ee_4,ee_3*ee_5}}

-- 1-dimensional subrepresentations
V1 =  ideal(w111 * matrix{{1},{1},{1}})
-- 2-dimensional subrepresentations
V2 =  ideal(w111 * matrix{{1,1},{xi,xi^2},{xi^2,xi}})

-- TEST: are these invariant?
assert isInvariant(V1)
assert isInvariant(V2)

-- the remaining monomials
EE2' = mingens ideal (EE2 % (V1+V2))

-- TEST: is this still invariant?
assert isInvariant(ideal EE2')

-- The subgroup S3 generated by D and E
sym3 = unique flatten apply(3,i->apply(3,j->D^i*E^j))

-- orbits on the remaining monomials
sym3orbits = unique apply(flatten entries EE2',ee->ideal mingens ideal apply(sym3,g->opE(ee,g)))

-- the dim 3 subrepresentations
V33 = select(sym3orbits,i->rank source mingens i==3)
V3 = V33#0
V3' = V33#1

-- TEST: are these invariant
assert isInvariant(V3)
assert isInvariant(V3')

-- the dim 6 subrepresentation
V6 = (select(sym3orbits,i->rank source mingens i==6))#0

-- TEST: are these invariant
assert isInvariant(V6)

-- the 9-dimensional representation from the paper
V = V2+V3+V6

-- the 3 dimensional representation from the paper
U = V1+V3'

-- the indices
ind = flatten apply(6,i->apply(i,j->{j,i}))

-- the Plucker space
R = K[apply(ind,i->x_i)]

-- the 2-forms in the order given by the index set
EE2R = matrix {apply(ind,i->ee_(i#0)*ee_(i#1))}

-- the representation on pluecker space
A2 = sub(opE(EE2R,A) // EE2R,K)
B2 = sub(opE(EE2R,B) // EE2R,K)
C2 = sub(opE(EE2R,C) // EE2R,K)
D2 = sub(opE(EE2R,D) // EE2R,K)
E2 = sub(opE(EE2R,E) // EE2R,K)
one2 = id_(K^15)

-- TEST: does this representation satisfy the relations?
assert (relationsG(A2,B2,C2,D2,E2,one2)==(14:true))

-- the operation on R
opR = (r,g) -> sub(r,(vars R)*g)

-- the skew matris
skew = matrix apply(6,i->apply(6,j->(
            if i<j then return x_{i,j};
            if i==j then return 0;
            if i>j then return -x_{j,i}
            )))

-- the plucker relations
betti res(G26 = ideal mingens pfaffians(4,skew))
-- 0: 1  .  .  .  .  . .
-- 1: . 15 35 21  .  . .
-- 2: .  .  . 21 35 15 .
-- 3: .  .  .  .  .  . 1

-- is the Grassmannian invariant
assert (opR(G26,A2) == G26)
assert (opR(G26,B2) == G26)
assert (opR(G26,C2) == G26)
assert (opR(G26,D2) == G26)
assert (opR(G26,E2) == G26)

-- the linear section
L = ideal ((vars R)*sub((gens U) // EE2R,K))

-- is the linear space invariant
assert (opR(L,A2) == L)
assert (opR(L,B2) == L)
assert (opR(L,C2) == L)
assert (opR(L,D2) == L)
assert (opR(L,E2) == L)

-- the section
betti (X = ideal mingens (G26+L))

-- TEST: is the variety invariant
assert (opR(X,A2) == X)
assert (opR(X,B2) == X)
assert (opR(X,C2) == X)
assert (opR(X,D2) == X)
assert (opR(X,E2) == X)


-- TEST: is the intersection transversal?
codim X  == codim G26 + codim L

-- the jacobina matrix
J = jacobian X

-- the coordinate points 
coordinatePoints = apply(15,i->matrix{apply(15,j->if i==j then 1 else 0)})

-- the coordinate points that lie on X
coordinatePointsOnX = select(coordinatePoints,point->0==sub(X,point))

-- The rank of the jacobian in these points
apply(coordinatePointsOnX,point-> rank sub(J,point))
-- {9, 9, 9, 9, 9, 9, 9, 9, 9}
-- all of these are singular points

-- find some points on X
randomPointsOnX = flatten apply(10,i->(
        betti (Xrandom = X + ideal(random(R^{0},R^{4:-1})));
        select(decompose Xrandom,i->1==degree i)
        )); #randomPointsOnX

-- TEST: are these indeed points?
apply(randomPointsOnX,i->assert (14==codim i));

-- is X smooth in these points?
apply(randomPointsOnX,i->rank sub(J,(vars R) % i))
-- {10, 10, 10, 10, 10, 10, 10, 10, 10, 10}

-- TEST: is X irreducible over K?
assert (1==#(decompose X))