(* ================================================================ *)

(*
Pfaffian2[L_] := (L[[1]][[2]] - L[[2]][[1]]) / 2
*)
Pfaffian2[L_] := L[[1]][[2]]

PfaffianHat[L_,j_] := 
	Drop[ Take[L,{2,Length[L]},{2,Length[L]}], {j-1}, {j-1} ]

(*
** Mathematica doesn't know if a is zero or not!!
** 1 == 0 is false; 0 == 0 is true; a == 0 evaluates to ... just a == 0.
PfaffianRecursive[L_] :=
	Sum[
		If[ L[[1]][[j]] == 1, 0,
			(-1)^j * L[[1]][[j]] * Pfaffian[PfaffianHat[L, j]]],
		{j, 2, Length[L]}]
*)

PfaffianRecursive[L_] :=
	Sum[(-1)^j * L[[1]][[j]] * Pfaffian[PfaffianHat[L, j]], {j, 2, Length[L]}]

(* xxx this currently doesn't check for even dims, nor for skew symmetry. *)
Pfaffian[L_] := If[ Length[L] == 2, Pfaffian2[L], PfaffianRecursive[L]]

(* ================================================================ *)
Submatrix2[L_,i_,j_] := {
	{ L[[i]][[i]], L[[i]][[j]] },
	{ L[[j]][[i]], L[[j]][[j]] }}

(* ================================================================ *)


PWLPW = {
{  0,  b,    a,  0,    0,  0,    0,  0 },
{ -b,  0,    0, -a,    0,  0,    0,  0 },
\
{ -a,  0,    0,  b,    a,  0,    0,  0 },
{  0,  a,   -b,  0,    0, -a,    0,  0 },
\
{  0,  0,   -a,  0,    0,  b,    a,  0 },
{  0,  0,    0,  a,   -b,  0,    0, -a },
\
{  0,  0,    0,  0,   -a,  0,    0,  b },
{  0,  0,    0,  0,    0,  a,   -b,  0 }
};

PWLPWinv = Simplify[Inverse[PWLPW]]

PbarLP = {
{  a,  0,    0,  0,    0,  0,    0,  0 },
{  0, -a,    0,  0,    0,  0,    0,  0 },
\
{  0,  0,    0,  0,    0,  0,   -a,  0 },
{  0,  0,    0,  0,    0,  0,    0,  a }
};

PWLPWdiff = {
{ -a,  0,     0,  0},
{  0,  a,     0,  0},
\
{  0,  0,     0,  0},
{  0,  0,     0,  0},
\
{  0,  0,     0,  0},
{  0,  0,     0,  0},
\
{  0,  0,     a,  0},
{  0,  0,     0, -a}
};

PWdiffLPW = {
{  a,  0,     0,  0,     0,  0,     0,  0 },
{  0, -a,     0,  0,     0,  0,     0,  0 },
\
{  0,  0,     0,  0,     0,  0,    -a,  0 },
{  0,  0,     0,  0,     0,  0,     0,  a }
};

X = Simplify[PWdiffLPW . PWLPWinv . PWLPWdiff];

Print["P_\\Omegadiff L P\\Omega"];      TeXForm[PWdiffLPW]
Print["(P_\\Omega L P_\\Omega)^{-1}"];  TeXForm[PWLPWinv]
Print["P_\\Omega  LP_\\Omegadiff"];     TeXForm[PWLPWdiff]
Print["X"]; TeXForm[X]

scale = Simplify[Pfaffian[-PWLPW]]
Print["Pf(-P_\\Omega LP_\\Omega)"]; TeXForm[scale]

(*
X={{ 0, a, b, c},
   {-a, 0, d, e},
   {-b,-d, 0, f},
   {-c,-e,-f, 0}};
*)  

Print["P_{11,00} X P_{11,00}="];  TeXForm[ Submatrix2[X,1,2]]
Print["Pf(-P_\\Omega L P_\\Omega) Pf(P_{11,00} X P_{11,00})="]
TeXForm[Simplify[scale * Pfaffian[ Submatrix2[X,1,2]]]]

Print["P_{10,10} X P_{10,10}="];  TeXForm[ Submatrix2[X,1,3]]
Print["Pf(-P_\\Omega L P_\\Omega) Pf(P_{10,10} X P_{10,10})="]
TeXForm[Simplify[scale * Pfaffian[ Submatrix2[X,1,3]]]]

Print["P_{10,01} X P_{10,01}="];  TeXForm[ Submatrix2[X,1,4]]
Print["Pf(-P_\\Omega L P_\\Omega) Pf(P_{10,01} X P_{10,01})="]
TeXForm[Simplify[scale * Pfaffian[ Submatrix2[X,1,4]]]]

Print["P_{01,10} X P_{01,10}="];  TeXForm[ Submatrix2[X,2,3]]
Print["Pf(-P_\\Omega L P_\\Omega) Pf(P_{01,10} X P_{01,10})="]
TeXForm[Simplify[scale * Pfaffian[ Submatrix2[X,2,3]]]]

Print["P_{01,01} X P_{01,01}="];  TeXForm[ Submatrix2[X,2,4]]
Print["Pf(-P_\\Omega L P_\\Omega) Pf(P_{01,01} X P_{01,01})="]
TeXForm[Simplify[scale * Pfaffian[ Submatrix2[X,2,4]]]]

Print["P_{00,11} X P_{00,11}="];  TeXForm[ Submatrix2[X,3,4]]
Print["Pf(-P_\\Omega L P_\\Omega) Pf(P_{00,11} X P_{00,11})="]
TeXForm[Simplify[scale * Pfaffian[ Submatrix2[X,3,4]]]]

Print["P_{11,11} X P_{11,11}= X"];
Print["Pf(-P_\\Omega L P_\\Omega) Pf(X)="]
TeXForm[Simplify[scale * Pfaffian[X]]]


(*
Simplify[Det[PWLPW ]]
Inverse[PWLPW ]*Det[PWLPW ]
TeXForm[Simplify[Det[PWLPW ]]]
TeXForm[Inverse[PWLPW ]*Det[PWLPW ]]
*)