-- computations for the article

-- "FORMATS OF 6×6 SKEW MATRICES OF LINEAR FORMS WITH VANISHING PFAFFIAN"

restart
R = QQ[m_{0,1}..m_{0,4},m_{1,2}..m_{1,4},a,b,c,y_0..y_2,l_0,l_1,l_2,l_{0,2}..l_{1,4}]

--------------------------
---- Proposition 1.4 -----
-- (computing pfaffians) -
--------------------------

Pf = (MM) -> (
    -- TEST: is MM skew?
    if not 0 == MM+transpose MM then error "Matrix not Skew";
    if rank source MM == 0 then return 1;
    -- special case i=0 
    i:=0;  
    -- since always i<= j
    theta:=0; 
    -- the laplacian analouge
    sum apply(rank source MM,j->(  
            -- the other indices
            notij = select(rank source MM,k->k != i and k !=j);         
            if i != j then (-1)^(i+j+theta+1)*MM_i_j*Pf(MM_notij^notij) else 0
            ))
    )

-- TEST: does this agree (up to sign) with the formula
-- Macaulay uses?
(
     Mtest0 = random(R^{6:0},R^{6:-1});
     Mtest = Mtest0 - transpose Mtest0;
     assert (ideal Pf(Mtest) == pfaffians(6,Mtest))
)     


---------------------
-- Proposition 2.4 --
--- (arrow format) --
---------------------

M = matrix{
      {       0, m_{0,1}, m_{0,2}, m_{0,3}, m_{0,4}, a},
      {-m_{0,1},       0, m_{1,2}, m_{1,3}, m_{1,4}, b},
      {-m_{0,2},-m_{1,2},       0,     y_0,     y_1, 0},
      {-m_{0,3},-m_{1,3},    -y_0,       0,     y_2, 0},
      {-m_{0,4},-m_{1,4},    -y_1,    -y_2,       0, 0},
      {      -a,      -b,       0,       0,       0, 0}
      }

-- TEST: is M skew?
assert (0==M+transpose M)


----------------
-- dim⟨a,b⟩=2 --
----------------

-- check formula for the Pfaffian
assert (Pf(M) == (b*Pf(M_{0,2,3,4}^{0,2,3,4}))-a*Pf(M_{1,2,3,4}^{1,2,3,4}))

----------------------------- 
-- dim⟨a,b⟩=2,c∈⟨y0,y1,y2⟩ --
-----------------------------

cc = sum apply(3,i->l_i*y_i)

S = id_(R^6) - matrix{
    {0,0,0,0,0,   0},
    {0,0,0,0,0,   0},
    {0,0,0,0,0, l_2},
    {0,0,0,0,0,-l_1},
    {0,0,0,0,0, l_0},
    {0,0,0,0,0,   0}
    }

N = S*M*transpose S    

-- TEST: is the formula for N as claimed?
assert (0 == 
    M - matrix{
    {     0,     0,a*l_2,-a*l_1,a*l_0,0},
    {     0,     0,b*l_2,-b*l_1,b*l_0,0},
    {-a*l_2,-b*l_2,    0,     0,    0,0},
    { a*l_1, b*l_1,    0,     0,    0,0},
    {-a*l_0,-b*l_0,    0,     0,    0,0},
    {     0,     0,    0,     0,    0,0}
    }
    -N
    )

-- the 4x4 submatrices
N05 = N_{1,2,3,4}^{1,2,3,4}
M05 = M_{1,2,3,4}^{1,2,3,4}

-- TEST: is the formula for Pf(N_05) as claimed?
assert (Pf(N05)==Pf(M05)-b*cc)

-- the other 4x5 submatrices
N15 = N_{0,2,3,4}^{0,2,3,4}
M15 = M_{0,2,3,4}^{0,2,3,4}

-- TEST: is it indeed similar for N_15?
assert (Pf(N15)==Pf(M15)-a*cc)

-------------------------------------------
-- dim⟨a, b⟩ = 2, c \not\in ⟨y0, y1, y2⟩ --
-------------------------------------------

M' = sub(M,flatten apply({0,1},i->apply({2,3,4},j->m_{i,j}=>l_{i,j}*c)))

-- TEST: is the formula for the Pfaffian of M' OK?
assert (
    Pf(M') == -c*(
    (l_{1,4}*y_0-l_{1,3}*y_1+l_{1,2}*y_2)*a-
    (l_{0,4}*y_0-l_{0,3}*y_1+l_{0,2}*y_2)*b
    )
    ) 


restart

R = QQ[l_0,l_3,l_4,m,n,a_0,a_3,a_4,aa_0,aa_3,aa_4,aaa_0,aaa_3,aaa_4]

--------------------------
---- Proposition 1.4 -----
-- (computing pfaffians) -
--------------------------

Pf = (MM) -> (
    -- TEST: is MM skew?
    if not 0 == MM+transpose MM then error "Matrix not Skew";
    if rank source MM == 0 then return 1;
    -- special case i=0 
    i:=0;  
    -- since always i<= j
    theta:=0; 
    -- the laplacian analouge
    sum apply(rank source MM,j->(  
            -- the other indices
            notij = select(rank source MM,k->k != i and k !=j);         
            if i != j then (-1)^(i+j+1)*MM_i_j*Pf(MM_notij^notij) else 0
            ))
    )

-- TEST: does this agree (up to sign) with the formula
-- Macaulay uses?
(
     Mtest0 = random(R^{6:0},R^{6:-1});
     Mtest = Mtest0 - transpose Mtest0;
     assert (ideal Pf(Mtest) == pfaffians(6,Mtest))
)     

---------------------
-- Proposition 3.5 --
---------------------



M = matrix{
    {   0,l_3,0,   0,l_0,m},
    {-l_3,  0,0,-l_0,  0,n},
    {   0,  0,0,  -m, -n,0},
    {   0,l_0,m,   0,l_4,0},
    {-l_0,  0,n,-l_4,  0,0},  
    {  -m, -n,0,   0,  0,0}
    }
    
-- TEST: is M skew?
assert (0==M+transpose M)

-- TEST: is Pf(M)=0?
assert (0==Pf(M))

-- the permutation
perm = {0,3,1,4,5,2} 

-- permute the matrix
Mperm = M_perm^perm

-- the symmetric 2x2 matrix
A = matrix{{l_3,l_0},{l_0,l_4}}

-- TEST: is A symmetric?
assert (0==A-transpose A)

-- the 2x2 identity matrix
id2 = id_(R^2)

-- the zero matrix
Z = 0*id2

-- TEST: does M have the claimed form?
assert (0==Mperm - (
    (       Z|       A|m*id2)||
    (      -A|       Z|n*id2)||
    ((-m*id2)|(-n*id2)|    0)
    )
    )

-- the other symmetric matrices
A' = sub(A,apply({0,3,4},i->l_i=>a_i))
A'' = sub(A,apply({0,3,4},i->l_i=>aa_i))
A''' = sub(A,apply({0,3,4},i->l_i=>aaa_i))

-- subsittution A = A'+ m A'' + n A'''
subA = apply({0,3,4},i->l_i=>a_i+m*aa_i+n*aaa_i)

-- TEST: is this the correct substitution
assert (sub(A,subA)==(A'+m*A''+n*A'''))

-- the base change matric
S = (
    (id2|  Z|  A''')||
    (  Z|id2|(-A'' ))||
    (  Z|  Z|  id2 )
    )

-- apply S
SMpermSt = S*Mperm*transpose(S)

-- use A = A'+mA''+nA'''
Mperm' = sub(SMpermSt,subA)

-- TEST: does this matrix have the claimed form?
assert (0==Mperm' - (
    (       Z|      A'|m*id2)||
    (     -A'|       Z|n*id2)||
    ((-m*id2)|(-n*id2)|    0)
    )
    )


-- a 2x2 matrix T
T = random(QQ^2,QQ^2)

-- the transpose
Tt = transpose T

-- the base change matrix
TTT = T++T++Tt^-1

-- TEST: Does T operate only on A'
TTT * Mperm' * transpose TTT == (
    (       Z| T*A'*Tt|m*id2)||
    (-T*A'*Tt|       Z|n*id2)||
    ((-m*id2)|(-n*id2)|    0)
    )

-- this checks the calculation only for a random T, not
-- for every T


------------------
-- special case --
----- | a b | ----
----- | b 0 | ----
------------------

subSpecial = {a_4=>0}
-- a_3 is a, a_0 is b

-- TEST: is A' of the claimed form
assert (sub(A',subSpecial) == matrix{{a_3,a_0},{a_0,0}})

-- the inverse permutation
perminv = {0,2,5,1,3,4}

-- TEST: is this indeed the inverse permutation
assert (perm_perminv == {0,1,2,3,4,5})
assert (perminv_perm == {0,1,2,3,4,5})

-- apply the imverse permutation
(sub(Mperm',subSpecial))_perminv^perminv

-- | 0    a_3 0 0    a_0 m |
-- | -a_3 0   0 -a_0 0   n |
-- | 0    0   0 -m   -n  0 |
-- | 0    a_0 m 0    0   0 |
-- | -a_0 0   n 0    0   0 |
-- | -m   -n  0 0    0   0 |

------------------
-- general case --
----- | a 0 | ----
----- | 0 b | ----
------------------

subGeneral = {a_0=>0}
-- a_3 is a, a_4 is b

-- TEST: is A' of the claimed form
assert (sub(A',subGeneral) == matrix{{a_3,0},{0,a_4}})

-- apply the permutation
MM' = (sub(Mperm',subGeneral))_perminv^perminv

-- rearrange the rows an columns again
perm2 = {0,1,5,2,3,4}

MM'_perm2^perm2
-- | 0    a_3 m 0 0    0   |
-- | -a_3 0   n 0 0    0   |
-- | -m   -n  0 0 0    0   |
-- | 0    0   0 0 -m   -n  |
-- | 0    0   0 m 0    a_4 |
-- | 0    0   0 n -a_4 0   |