// Approaching a definitive form
// CONTENTS. The different forms of the syzygy matrix for
// FUN, Section 4. In order, they are
// p.25, Mz1, My1, Mx1 and Mla.
// All of this works easily in the online Magma calculator
// magma.maths.usyd.edu.au/calc

RRgr<la,mu,a,b,c,d,e,x1,x2,y1,y2,z1,z2> :=
PolynomialRing(Rationals(),[0,0,2,2,1,1,4,1,1,2,2,3,3]);
J16 := ZeroMatrix(RRgr,16,16); for i in [1..8] do J16[i,i+8] := 1;
  J16[i+8,i] := 1; end for;  // standard quadratic form 0 I \\ I 0

// the form of the 9 equations given on p. 25 of Fun2
M := Matrix(9,[
y1, -x1, d*x2, la, -mu*d, 0, 0, 0, 0,
c*x1, -x2, y2, 0, la*c, mu, 0, 0, 0,
-mu*b*d, la*b, z1, -x1, y1, 0, mu, 0, 0,
la*b*c, 0, z2, -x2, c*x1, 0, 0, mu, 0,
0, z1, -la*e, -d*x2, 0, -y1, la*c, mu*d, 0,
-mu*e, z2, 0, -y2, 0, -c*x1, 0, la*c, mu*d,
z1, 0, mu*a*d, 0, -d*x2, -x1, 0, la, 0,
z2, mu*a, -la*a*c, 0, -y2, -x2, 0, 0, la,
0, la*a*b*c, mu*a*e, -a*c*x1, -e*x2, b*y2, 0, -z2, z1,
la*b*e, mu*a*b*d, 0, a*y1, e*x1, -b*d*x2, z2, -z1, 0,
e*x2, a*c*x1, 0, -la*a*c, 0, -z2, 0, y2, -d*x2,
-e*x1, -a*y1, 0, -mu*a*d, -la*e, z1, -y2, d*x2, 0,
b*y2, 0, -a*c*x1, mu*a, -z2, 0, 0, -x2, x1,
-b*d*x2, 0, a*y1, 0, z1, la*b, x2, -x1, 0,
0, -b*y2, -e*x2, z2, mu*e, -la*b*c, 0, c*x1, -y1,
0, b*d*x2, e*x1, -z1, 0, -mu*b*d, -c*x1, y1, 0
]);

L := [
    -mu^2*a*d + d*x2^2 - x1*y2 + la*z2,
    la*mu*e + c*d*x1*x2 - y1*y2 + la*c*z1 + mu*d*z2,
    -la^2*b*c + c*x1^2 - x2*y1 + mu*z1,
    mu*e*x1 + la*b*c*d*x2 + mu*b*d*y2 + c*x1*z1 - y1*z2,
    -mu*a*y1 + la*b*y2 + x2*z1 - x1*z2,
    mu*a*c*d*x1 + la*e*x2 + la*a*c*y1 - y2*z1 + d*x2*z2,
    -mu^2*a*b*d^2 - la^2*b*e + e*x1^2 + b*d^2*x2^2 + a*y1^2 - z1^2,
    la*mu*a*b*c*d + e*x1*x2 + a*c*x1*y1 + b*d*x2*y2 - z1*z2,
    -la^2*a*b*c^2 - mu^2*a*e + a*c^2*x1^2 + e*x2^2 + b*y2^2 - z2^2
];

// Basic sanity checks
Matrix(9,L)*Transpose(M);
Transpose(M)*J16*M;

// Calculate some of the spinors
// Determinant(Submatrix(M,[1,2,3,4,5,6,7,16],[1..8])) div L[9];
// SquareRoot(-$1) div (mu*d); L[3];

for i in [1..8] do
 si := SquareRoot(-(-1)^i*Determinant(Submatrix(M, Remove([1..8],i) cat [8+i], [1..8])) div L[9]);
si;
Factorization(si);
end for;


// Matrix Mz1
Mz1 := Matrix(9, [
0, a*c*x1, b*y2, z2, 0, mu*a, x1, x2, 0,
-a*c*x1, 0, x2*e, 0, z2, -la*a*c, -d*x2, -y2, 0,
-b*y2, -x2*e, 0, -mu*e, la*b*c, z2, y1, c*x1, 0,
-z2, 0, mu*e, 0, -c*x1, -y2, -mu*d, la*c, 0,
0, -z2, -la*b*c, c*x1, 0, -x2, 0, mu, 0,
-mu*a, la*a*c, -z2, y2, x2, 0, la, 0, 0,
-x1, d*x2, -y1, mu*d, 0, -la, 0, 0, 0,
-x2, y2, -c*x1, -la*c, -mu, 0, 0, 0, 0,
z1, la*e, 0, 0, y1, d*x2, 0, -mu*d, la*c,
-la*b, z1, -mu*b*d, -y1, 0, x1, 0, 0, mu,
0, mu*a*d, z1, -d*x2, -x1, 0, 0, la, 0,
0, a*y1, b*d*x2, z1, -la*b, 0, 0, -x1, -x2,
-a*y1, 0, x1*e, la*e, z1, mu*a*d, 0, d*x2, y2,
-b*d*x2, -x1*e, 0, 0, -mu*b*d, z1, 0, -y1, -c*x1,
la*a*b*c, mu*a*e, 0, -x2*e, b*y2, -a*c*x1, z1, z2, 0,
mu*a*b*d, 0, -la*b*e, x1*e, -b*d*x2, a*y1, 0, z1, z2
]);

Lz1 := [c*d*x1*x2-y1*y2+la*c*z1+mu*d*z2-la*mu*e,
la^2*b*c+c*x1^2-x2*y1+mu*z1,
-mu^2*a*d-d*x2^2+x1*y2-la*z2,
-mu*a*y1+la*b*y2-x2*z1+x1*z2,
mu*a*c*d*x1+la*a*c*y1+y2*z1-d*x2*z2+la*x2*e,
la*b*c*d*x2+mu*b*d*y2-c*x1*z1+y1*z2+mu*x1*e,
-la^2*a*b*c^2-a*c^2*x1^2-b*y2^2-z2^2-mu^2*a*e-x2^2*e,
-la*mu*a*b*c*d+a*c*x1*y1+b*d*x2*y2+z1*z2+x1*x2*e,
-mu^2*a*b*d^2-b*d^2*x2^2-a*y1^2-z1^2-la^2*b*e-x1^2*e];

// sanity checks
Matrix(9,Lz1)*Transpose(Mz1);
Transpose(Mz1)*J16*Mz1;

// i -> 3 and i -> 7 give relations for x1y2 and z2^2 as common factors.
i:=3; [Reversion(Factorization(m))[1,1] :
 m in  Pfaffians(Submatrix(Mz1,I,I),6) where I is Remove([1..8],i) | m ne 0];

// typical SO(16) row operation
AddRow(AddRow(Mat,-a*c,7,1),a*c,9,15);
// that is, add X times Row i to Row j and -X times Row j+8 to row i+8 (mod 16).

> Pfaffian(Submatrix(Mz1,[1..6],[1..6])) div z2;
// -la^2*a*b*c^2 - a*c^2*x1^2 - b*y2^2 - z2^2 - mu^2*a*e - x2^2*e

> Factorization(Determinant(Submatrix(Mz1,[9..14],[1..6])));
// z1^2*(mu^2*a*b*d^2 + b*d^2*x2^2 + a*y1^2 + z1^2 + la^2*b*e + x1^2*e)^2

for i in [1..7], j in [i+1..8] do [i,j]; Reversion(Factorization(Determinant(Submatrix(Mz1,J,I))))[1] where J is [l+8 : l in I] where I is Remove(Remove([1..8],j),i); end for;
for i in [1..7], j in [i+1..8] do [i,j]; I := Remove(Remove([1..8],j),i); X := Pfaffian(Submatrix(Mz1,I,I)); if X eq 0 then 0; else Factorization(X); end if; end for;


// Matrix My1

My1 := Matrix(9,[
0, a*c*x1, x2*e, b*y2, mu*a*e, z2, z1, la*a*b*c, 0,
-c*x1, 0, la*c, -mu, y2, 0, 0, -x2, 0,
-x2*e, -la*a*c, 0, -z2, 0, y2, d*x2, a*c*x1, 0,
-b*y2, mu*a, z2, 0, -a*c*x1, -x2, -x1, 0, 0,
-mu*e, -y2, 0, c*x1, 0, -la*c, mu*d, z2, 0,
-z2, 0, -y2, x2, la*a*c, 0, la, -mu*a, 0,
-z1, 0, -d*x2, x1, -mu*a*d, -la, 0, 0, 0,
-la*b*c, x2, -c*x1, 0, -z2, mu, 0, 0, 0,
-y1, la, -mu*d, 0, d*x2, 0, 0, -x1, 0,
-la*b*e, -a*y1, -x1*e, -b*d*x2, 0, z1, 0, mu*a*b*d, z2,
mu*b*d, x1, -y1, 0, -z1, 0, 0, la*b, -mu,
0, d*x2, 0, -y1, -la*e, mu*d, 0, -z1, -la*c,
-b*d*x2, 0, z1, la*b, -a*y1, x1, 0, 0, x2,
0, -z1, 0, -mu*b*d, -x1*e, -y1, 0, -b*d*x2, -c*x1,
0, z2, -mu*e, -la*b*c, x2*e, -c*x1, -y1, b*y2, 0,
x1*e, -mu*a*d, -la*e, z1, 0, d*x2, 0, -a*y1, y2
]);

My1a := MultiplyRow(MultiplyRow(MultiplyRow(M13,a,2),a,5),a,8);
Pfaffians(Submatrix(My1a,I,I),6) where I is [1,2,3,4,5,6,8];
[Reversion(Factorization(m))[1,1] : m in $1 | m ne 0];

Factorization(Determinant(Submatrix(My1a,J,I))) where J is [i+8 : i in I] where I is Remove([1..8],7);
i := 1; [Factorization(m) : m in Pfaffians(Submatrix(My1a,I,I),6) | m ne 0] where I is Remove([1..8],i);


Mx1 := Matrix(9,[
0, z1, y1, d*x2, mu*d, la*c, la*b*c*d, 0, la*e,
-z1, 0, -b*y2, 0, z2, -mu*a*c, 0, a*c*y1+x2*e, -mu*a*e,
-y1, b*y2, 0, -z2, 0, c*x2, 0, -c*z1+mu*e, x2*e,
-d*x2, 0, z2, 0, y2, 0, a*c*y1+x2*e, -mu*a*c*d, 0,
-mu*d, -z2, 0, -y2, 0, 0, -c*z1+mu*e, c*d*x2, 0,
-la, mu*a, -x2, 0, 0, 0, z2, -y2, 0,
-la*b*d, 0, 0, -a*y1, z1, -z2, la*b*e, -b*d*y2-x1*e, 0,
0, -a*y1, z1, mu*a*d, -d*x2, y2, b*d*y2+x1*e, la*e, 0,
x1, 0, 0, mu*a, -x2, 0, b*y2, z2, 0,
0, x1, 0, la, 0, 0, y1, -mu*d, -d*x2,
0, 0, x1, 0, la, 0, -z1, d*x2, -mu*a*d,
0, -la*b, 0, x1, 0, mu, 0, -y1, z1,
0, 0, -la*b, 0, x1, -x2, 0, z1, a*y1,
0, 0, 0, -z1, -y1, c*x1, 0, -la*b*c*d, b*d*y2+x1*e,
0, x2, mu, 0, 0, 0, c*x1, la*c, -y2,
0, 0, 0, -x2, -mu, 0, -la*b*c, c*x1, -z2
]);

Lx1 := [a*c*x2*y1+b*y2^2-mu*a*c*z1+z2^2+mu^2*a*e+x2^2*e,
-c*d*x1*x2+y1*y2-la*c*z1-mu*d*z2+la*mu*e,
-mu*a*c*d*x1-la*a*c*y1-y2*z1+d*x2*z2-la*x2*e,
-la*b*c*d*x2-mu*b*d*y2+c*x1*z1-y1*z2-mu*x1*e,
-la*mu*a*b*c*d+a*c*x1*y1+b*d*x2*y2+z1*z2+x1*x2*e,
a*y1^2+b*d*x1*y2+z1^2-la*b*d*z2+la^2*b*e+x1^2*e,
mu^2*a*d+d*x2^2-x1*y2+la*z2,
mu*a*y1-la*b*y2+x2*z1-x1*z2,
-la^2*b*c-c*x1^2+x2*y1-mu*z1];


Mla := Matrix(9,[
0, -y1, -z1, c*x1, -b*d*x2, -mu*b*d, 0, b*c*d*x1, x1*e,
y1, 0, -z2, c*x2, -b*y2, 0, -c*z1+mu*e, 0, x2*e,
z1, z2, 0, -mu*a*c, 0, b*y2, a*c*y1+x2*e, 0, -mu*a*e,
-x1, -x2, mu*a, 0, 0, 0, z2, b*y2, 0,
d*x2, y2, 0, 0, 0, -z2, -mu*a*c*d, a*c*y1+x2*e, 0,
mu*d, 0, -y2, 0, z2, 0, c*d*x2, -c*z1+mu*e, 0,
0, z1, -a*y1, -z2, mu*a*b*d, -b*d*x2, -x1*e, -b*d*z2+la*b*e, 0,
-d*x1, 0, 0, -y2, -a*y1, z1, d*z2-la*e, -x1*e, 0,
la, 0, 0, 0, -mu*a, x2, -y2, z2, 0,
0, la, 0, 0, 0, -x1, d*x2, -z1, -mu*a*d,
0, 0, la, 0, -x1, 0, -mu*d, y1, -d*x2,
0, 0, 0, la*c, -z1, -y1, c*d*x1, 0, -d*z2+la*e,
0, 0, x1, mu, la*b, 0, -y1, 0, z1,
0, x1, 0, -x2, 0, la*b, z1, 0, a*y1,
0, 0, 0, 0, -x2, -mu, la*c, c*x1, -y2,
0, mu, x2, 0, 0, 0, -c*x1, la*b*c, z2
]);

Lla :=
[a*c*x2*y1 + b*y2^2 - mu*a*c*z1 + z2^2 + mu^2*a*e + x2^2*e,
la*mu*a*b*c*d - a*c*x1*y1 - b*d*x2*y2 - z1*z2 - x1*x2*e,
la*b*c*d*x2 + mu*b*d*y2 - c*x1*z1 + y1*z2 + mu*x1*e,
-a*y1^2 - b*d*x1*y2 - z1^2 + la*b*d*z2 - la^2*b*e - x1^2*e,
-c*d*x1*x2 + y1*y2 - la*c*z1 - mu*d*z2 + la*mu*e,
-mu*a*c*d*x1 - la*a*c*y1 - y2*z1 + d*x2*z2 - la*x2*e,
-mu*a*y1 + la*b*y2 - x2*z1 + x1*z2,
-mu^2*a*d - d*x2^2 + x1*y2 - la*z2,
la^2*b*c + c*x1^2 - x2*y1 + mu*z1];

// to make its top left block skew:
> AddRow(AddRow(MultiplyRow(MultiplyRow(MultiplyRow(MultiplyRow(MultiplyRow(M1,c,4),b,5),b,6),c,7),b*c,8),-e,16,7),b*e,15,8);
> M + Transpose(M) where M is Submatrix($1,[1..8],[1..8]);