//////////////////////////////////////////////////////////////////////
// Computing the cascading Pfaffians of a flip
//////////////////////////////////////////////////////////////////////

/*
How to use
==========
Find the file 'lr.m' in directory PATH and attach it to magma during
a magma session  by doing

&gt; Attach("PATH/lr.m");

Compute the pfaffian equations of a Type I long rectangle pair
as follows.  The input must be a sequence of the form [d,e,d,e,...]
(of any length) with d,e &gt;= 2.

&gt; X := TypeIPair([4,3,4]);
&gt; Print(X);
        4       1       0       -2
        ----------      ----------
        3       2       3       2
        4       2       4       2
                3               3
                2               2
        ----------      ----------
        0       -3      3       1

&gt; E := PfaffianEquations(X);
&gt; E;
{
    -x1^3*A^2*B^7 - y0*M + x0*y1,
    -A*B^4*M - x1^4 + x0*x2,
    -x2^2*A*B^3 - x0*L*M^2 + x1*y2,
    -B*L*M^3 - x2^3 + x1*x3,
    -x2*A^2*B^7 - x0^2*L*M + x1*y1,
    -x1^2*L^2*M^5 - x3^2*A*B^2 + x2*y3,
    -y2*L*M^3 - x2^2*x3*A*B^2 + x1*y3,
    -x1*L^3*M^8 - x3^3*A*B + x2*y4,
    -x3*A*B^3 - x1^3*L*M^2 + x2*y2,
    -y3*L*M^3 - x2^2*x3^2*A*B + x1*y4,
    -x1^3*x2*A*B^3 - y1*M + x0*y2,
    -y4*L*M^3 - x2^2*x3^3*A + x1*y5,
    -y1*A*B^4 - x0^2*x1^3*L + x2*y0,
    -A^3*B^11 - x0^3*L + x1*y0,
    -y2*B*M - x1^3*x2^2 + x0*x3,
    -L^4*M^11 - x3^4*A + x2*y5,
    -x2^2*L^3*M^8 - y5*B + x3*y4
}

The equations are not written in any particular order.  But the
terms of each are weighted so that (with luck) they are written
in the reverse order to that done by hand.

See the working of this calculation by doing

&gt; SetVerbose("User1",true);   
&gt; E := PfaffianEquations(X);

STEP 1:
        square: x3,x2,y5,y4 with equations
                x2*y5 = x3^4*A + L^4*M^11
                y4*x3 = y5*B + x2^2*L^3*M^8
STEP 2: adjoin x1
        pentagram x1,x2,y4,y5,x3 with matrix
[      0      x2   L*M^3  x3^3*A      x1]
[      0       0      x3 L^3*M^8       B]
[      0       0       0      y5    x2^2]
[      0       0       0       0      y4]
        eliminate y5
        square x1,y4,x3,x2 with equations
                x1*x3 = x2^3 + B*L*M^3
                y4*x2 = x3^3*A*B + x1*L^3*M^8
STEP 3: adjoin y3
        pentagram y3,x1,y4,x3,x2 with matrix
[         0         x1          B       x2^2         y3]
[         0          0         x2      L*M^3   x3^2*A*B]
[         0          0          0         x3 x1*L^2*M^5]
[         0          0          0          0         y4]
        eliminate y4
        square x1,y3,x3,x2 with equations
                x1*x3 = x2^3 + B*L*M^3
                y3*x2 = x3^2*A*B^2 + x1^2*L^2*M^5
STEP 4: adjoin y2
        pentagram y2,x1,y3,x3,x2 with matrix
[         0         x1          B       x2^2         y2]
[         0          0         x2      L*M^3   x3*A*B^2]
[         0          0          0         x3 x1^2*L*M^2]
[         0          0          0          0         y3]
        eliminate y3
        square x1,y2,x3,x2 with equations
                x1*x3 = x2^3 + B*L*M^3
                y2*x2 = x3*A*B^3 + x1^3*L*M^2
STEP 5: adjoin x0
        pentagram x0,x1,y2,x3,x2 with matrix
[    0    x1   B*M  x2^2    x0]
[    0     0    x2 L*M^2 A*B^3]
[    0     0     0    x3  x1^3]
[    0     0     0     0    y2]
        eliminate x3
        square x0,y2,x2,x1 with equations
                x0*x2 = x1^4 + A*B^4*M
                y2*x1 = x2^2*A*B^3 + x0*L*M^2
STEP 6: adjoin y1
        pentagram y1,x0,y2,x2,x1 with matrix
[       0       x0    A*B^4     x1^3       y1]
[       0        0       x1        M x2*A*B^3]
[       0        0        0       x2   x0*L*M]
[       0        0        0        0       y2]
        eliminate y2
        square x0,y1,x2,x1 with equations
                x0*x2 = x1^4 + A*B^4*M
                y1*x1 = x2*A^2*B^7 + x0^2*L*M
STEP 7: adjoin y0
        pentagram y0,x0,y1,x2,x1 with matrix
[      0      x0   A*B^4    x1^3      y0]
[      0       0      x1       M A^2*B^7]
[      0       0       0      x2  x0^2*L]
[      0       0       0       0      y1]
        eliminate x2
        square y0,x0,x1,y1 with equations
                y0*x1 = A^3*B^11 + x0^3*L
                x0*y1 = x1^3*A^2*B^7 + y0*M

Switch the printing off by

&gt; SetVerbose("User1",false);   

*/


//////////////////////////////////////////////////////////////////////
//		Defining a pair of long rectangles
//////////////////////////////////////////////////////////////////////

LR := recformat&lt;		name	: MonStgElt,
	left, right, proj, zerocf	: SeqEnum,
				top 	: BoolElt	&gt;; 

LRP := recformat&lt;		name	: MonStgElt,
		 	left, right	: Rec,
			 	ring	: RngMPol	&gt;; 

intrinsic TypeIPair(L::SeqEnum) -&gt; Rec
{A pair of long rectangles with short continued fraction [d,e,d,...]
with d,e &gt;= 2.}
	require #L ge 2 and &amp;and[ d ge 2 : d in L ]:
		"Illegal argument; must be [d,e,d,...] with d,e &gt;= 2";
	for i in [3..#L by 2] do
		require L[i] eq L[1]:
		"Illegal argument; must be [d,e,d,...] with d,e &gt;= 2";
	end for;
	for i in [4..#L by 2] do
		require L[i] eq L[2]:
		"Illegal argument; must be [d,e,d,...] with d,e &gt;= 2";
	end for;
	Xl := TypeILongRectangle(L);
	if L[1] eq L[#L] then
		R := L[2..#L] cat [L[2]];
	else
		R := L;
	end if;
	Xr := TypeILongRectangle(R : Top:=false);
	nx := #Xr`left;
	ny := #Xr`right;
	n := nx + ny + 4;
	M := IdentityMatrix(Rationals(),n);
	firstrows := Eltseq(M)[1..n*(n-2)];
	monomialorderweights :=
		[1:i in [1..n-4]] cat [100,100,100,100]
		cat [ i : i in [1..nx] ] cat [ i : i in [1..ny] ]
				cat [100,200,300,400]
		cat firstrows;
	ring := PolynomialRing(Rationals(),n,"weight",monomialorderweights);
	X := rec&lt; LRP |	name := "pair",
		 left := Xl, right := Xr, ring := ring &gt;;
		// name the variables on X
	names := [];
	for i in [1..nx] do
		Append(~names,"x" cat IntegerToString(i-1));
	end for;
	for i in [1..ny] do
		Append(~names,"y" cat IntegerToString(i-1));
	end for;
	names cat:= ["A", "B", "L", "M"];
	AssignNames(~ring,names);
	return X;
end intrinsic;

intrinsic TypeILongRectangle(L::SeqEnum : Top:=true) -&gt; Rec
{}
	if Top then
		left := L cat [0];
		ccf := ComplementaryContinuedFraction(L[1..#L-1]);
		right := [1] cat ccf cat [-L[#L]+1];
		zerocf := Reversion(L[1..#L-1]) cat [1] cat ccf;
	else
		left :=[0] cat Reversion(L);
		ccf := ComplementaryContinuedFraction(L[1..#L-1]);
		right := [-L[#L]+1] cat Reversion(ccf) cat [1];
		zerocf := Reversion(L[1..#L-1]) cat [1] cat ccf;
	end if;
	X := rec&lt; LR | left:=left, right:=right, zerocf:=zerocf,
		top := Top	&gt;;
	return X;
end intrinsic;

intrinsic Left(X::Rec) -&gt; Rec
{}
	return X`left;
end intrinsic;

intrinsic Right(X::Rec) -&gt; Rec
{}
	return X`right;
end intrinsic;


//////////////////////////////////////////////////////////////////////
// Printing
//////////////////////////////////////////////////////////////////////

intrinsic Print(X::Rec)
{}
	require assigned X`name:
		"Argument is not a long rectangle (pair)";
	name := X`name;
	require name in {"pair","lr"}:
		"Argument is not a long rectangle (pair)";
	if name eq "pair" then
		Xl := X`left;
		Xr := X`right;
		Ll := Xl`left;
		Rl := Xl`right;
		Lr := Xr`left;
		Rr := Xr`right;
		L := Ll;
		R := Rl;
		printf "\t%o\t%o",Ll[1],Rl[1];
		printf "\t%o\t%o\n",Lr[1],Rr[1];
		printf "\t----------\t----------\n";
		for i in [2..Maximum(#L,#R)] do
			if i le #L-1 and i le #R-1 then
				printf "\t%o\t%o",Ll[i],Rl[i];
				printf "\t%o\t%o\n",Lr[i],Rr[i];
			elif i gt #L-1 and i le #R-1 then
				printf "\t\t%o",Rl[i];
				printf "\t\t%o\n",Rr[i];
			elif i le #L-1 and i gt #R-1 then
				printf "\t%o\t",Ll[i];
				printf "\t%o\t\n",Lr[i];
			end if;
		end for;
		printf "\t----------\t----------\n";
		printf "\t%o\t%o\t%o\t%o\n",Ll[#L],Rl[#R],Lr[#L],Rr[#R];
	elif name eq "lr" then
		L := X`left;
		R := X`right;
		printf "\t%o\t%o\n\t----------\n",L[1],R[1];
		for i in [2..Maximum(#L,#R)] do
			if i le #L-1 and i le #R-1 then
				printf "\t%o\t%o\n",L[i],R[i];
			elif i gt #L-1 and i le #R-1 then
				printf "\t\t%o\n",R[i];
			elif i le #L-1 and i gt #R-1 then
				printf "\t%o\t\n",L[i];
			end if;
		end for;
		printf "\t----------\n\t%o\t%o\n",L[#L],R[#R];
	end if;
end intrinsic;


//////////////////////////////////////////////////////////////////////
//		Computing the cascading Pfaffians of a pair
//////////////////////////////////////////////////////////////////////

forward makepfaffian, is_in, is_bottom, addeqn, projections;
intrinsic PfaffianEquations(X::Rec) -&gt; SeqEnum
{A sequence containing the (encoded) Pfaffian equations of X, a pair
of long rectangles}
	Pl := projections(X`left);	// records top vars and annotations
	Pr := projections(X`right);	// records bottom vars and annotations
	R := X`ring;			// poly ring with all vars
	ABleft := Left(X)`left;		// seqs ctg the side tags
	ABright := Left(X)`right;
	LMleft := Right(X)`left;
	LMright := Right(X)`right;
	nx := #ABleft;
	ny := #ABright;

		//	*** FIRST STEP ***
	// Equations are of the form
	//	x * y = term1 + term2,  where L,M appear only in term2.
	// This is encoded in a sequence as
	//	[ x, y, term1, term2 ]
	// where x,y are ordered so that x is at the bottom end of
	// the four variables concerned.
	// We also order the equations so that E1 starts with the
	// lowest left-hand variable.
	x1 := R.(Pl[1][1]+1);
	x0 := R.(Pr[1][1]+1);
	y1 := R.(Pl[1][2]+nx+1);
	y0 := R.(Pr[1][2]+nx+1);
	n := nx+ny;
	A := R.(n+1);
	B := R.(n+2);
	L := R.(n+3);
	M := R.(n+4);
	equations := { R | };
	pl := Pl[1];		// eg [3,4,1,1,2,1] means x1^3AB, y1^4A^2B
	pr := Pr[1];		// eg [1,2,1,0,0,1] means x0L, y0^2M
	E1 := [x0,y1, x1^ABleft[1]*A^pl[3]*B^pl[4],
		y0^(-LMleft[1])*L^pr[5]*M^pr[6]];
	E2 := [y0,x1, y1^ABright[1]*A^pl[5]*B^pl[6],
		x0^(-LMright[1])*L^pr[3]*M^pr[4]];
	addeqn(~equations,E1);
	addeqn(~equations,E2);
	vprintf User1: "STEP 1:\n";
	vprintf User1: "\tsquare: %o,%o,%o,%o with equations\n",x1,x0,y1,y0;
	vprintf User1: "\t\t%o*%o = %o + %o\n",E1[1],E1[2],E1[3],E1[4];
	vprintf User1: "\t\t%o*%o = %o + %o\n",E2[1],E2[2],E2[3],E2[4];

		//	*** MAIN LOOP ***
	for i in [2..#Pl] do
		// work out new variable at the bottom, then make the
		// pfaffian with new var and two previous equations
		plold := pl;
		prold := pr;
		pl := Pl[i];
		pr := Pr[i];
		vprintf User1: "STEP %o:\t",i;
		if pr[1] eq prold[1] then		// this finds new var
			newvar := R.(pr[2]+nx+1);
		else
			newvar := R.(pr[1]+1);
		end if;
		// new pentagram is easy to write down now
		penta := [ newvar, E1[1],E2[1],E1[2],E2[2] ];
		vprintf User1:
			 "adjoin %o\n\tpentagram %o,%o,%o,%o,%o with matrix\n",
			newvar,penta[1],penta[2],penta[3],penta[4],penta[5];
		M := makepfaffian(E1,E2,newvar);	// builds the matrix
		vprintf User1: "%o\n",M;
		// compute three new equations in the sequence format
		// (but don't worry about the x,y order yet)
		E3 := [ M[1,5], M[2,3], M[1,2]*M[3,5], M[1,3]*M[2,5] ];
		E4 := [ M[1,2], M[4,5], M[1,4]*M[2,5], M[1,5]*M[2,4] ];
		E5 := [ M[1,5], M[3,4], M[1,4]*M[3,5], M[1,3]*M[4,5] ];

		// eliminate var from top to leave two equations
		// making up the new 4-var square
		if pl[1] eq plold[1] then
			elimvar := R.(plold[2]+nx+1);
		else
			elimvar := R.(plold[1]+1);
		end if;
		// find and normalise the next two equations;  they are
		// the 2 of 5 we own that do not involve elimvar
		vprintf User1: "\teliminate %o",elimvar;
		eqns := [ E1, E2, E3, E4, E5 ];
		neweqns := [];
		for E in eqns do
			if not is_in(elimvar,E) then
				Append(~neweqns,E);
			end if;
		end for;
		E1 := neweqns[1];
		E2 := neweqns[2];
		Espare := [ E : E in [E3, E4, E5] | E ne E1 and E ne E2 ][1];
		/* this next bit is wrong:  sometimes a Pfaffian eqn
			does not appear if we do it this way
		addeqn(~equations,Espare);
		addeqn(~equations,E1);
		addeqn(~equations,E2);
		*/
		// as overkill, record all the pfaffian equations
		addeqn(~equations,eqns[1]);
		addeqn(~equations,eqns[2]);
		addeqn(~equations,eqns[3]);
		addeqn(~equations,eqns[4]);
		addeqn(~equations,eqns[5]);
		// finally, sort the first two entries of E1,E2 ...
		fourvars := &amp;cat[ E1[1..2], E2[1..2] ];
		if not is_bottom(E1[1],penta,nx) then
			Etemp := E1;
			E1[1] := Etemp[2];
			E1[2] := Etemp[1];
		end if;
		if not is_bottom(E2[1],penta,nx) then
			Etemp := E2;
			E2[1] := Etemp[2];
			E2[2] := Etemp[1];
		end if;
		// ... and sort the last two entries:  L or M must divide
		//  the 4th entry
		if IsDivisibleBy(E1[3],L) or IsDivisibleBy(E1[3],L) then
			Etemp := E1;
			E1[3] := Etemp[4];
			E1[4] := Etemp[3];
		end if;
		if IsDivisibleBy(E2[3],L) or IsDivisibleBy(E2[3],L) then
			Etemp := E2;
			E2[3] := Etemp[4];
			E2[4] := Etemp[3];
		end if;
		vprintf User1:
			"\n\tto get\t%o*%o = %o + %o",
				Espare[1],Espare[2],Espare[3],Espare[4];
		vprintf User1:
			"\n\tand the square %o,%o,%o,%o with equations\n",
				E1[1],E2[1],E1[2],E2[2];
		vprintf User1: "\t\t%o*%o = %o + %o\n",E1[1],E1[2],E1[3],E1[4];
		vprintf User1: "\t\t%o*%o = %o + %o\n",E2[1],E2[2],E2[3],E2[4];
	end for;
	return equations;
end intrinsic;


//////////////////////////////////////////////////////////////////////
//			Auxilliary functions
//////////////////////////////////////////////////////////////////////

intrinsic ComplementaryContinuedFraction(Q::[RngIntElt]) -&gt; SeqEnum
{}
	require &amp;and [ q ge 2 : q in Q ]: "Need all numbers &gt;= 2";
	result := [ Integers() | ];
	for i in [1..&amp;+Q] do
		if Q[1] eq 1 then
			if #Q eq 1 then break; end if;
			Q := Q[2..#Q];
			m := result[#result];
			m +:= 1;
			result[#result] := m;
		else
			Append(~result,2);
		end if;
		n := Q[1];
		n -:= 1;
		Q[1] := n;
	end for;
	return result;
end intrinsic;

procedure addeqn(~E,f)
nn := #E;
	f1 := f[1]*f[2] - f[3] - f[4];
	Include(~E,f1);
end procedure;

Rindex := func&lt; R,x | Index([R.i:i in [1..Rank(R)]],x) &gt;;

is_x := func&lt; x, nx | Rindex(Parent(x),x) le nx &gt;;

function is_bottom(x,penta,nx)
	R := Parent(x);
	xinds := [ Rindex(R,y) : y in penta | is_x(y,nx) ];
	yinds := [ Rindex(R,y) : y in penta | not is_x(y,nx) ];
	if is_x(x,nx) then
		return Rindex(R,x) eq Minimum(xinds);
	else
		return Rindex(R,x) eq Minimum(yinds);
	end if;
end function;

function is_in(x,E)
	return &amp;or [ IsDivisibleBy(E[i],x) : i in [1..4] ];
end function;

function makepfaffian(E1,E2,x)
	a15 := x;
	a12 := E1[1];
	a34 := E1[2];
	a23 := E2[2];
	a45 := E2[1];
	a25 := E2[3] div a34;
	a14 := E1[3] div a23;
	a24 := GCD(E1[4],E2[4]);
	a13 := E1[4] div a24;
	a35 := E2[4] div a24;
	M := Matrix([	[0,a12,a13,a14,a15],
			[0,0,a23,a24,a25],
			[0,0,0,a34,a35],
			[0,0,0,0,a45]]);
	return M;
end function;

function projections(X)
	if assigned X`proj then
		return X`proj;
	end if;
	L := X`left;
	R := X`right;
	if not X`top then
		result := [ [0,0,1,0,0,1] ];
		while #L gt 2 or #R gt 2 do
			ij := result[#result];
			if (L[#L] eq 1 and R[#R] eq 1 and #L eq 3) or
							L[#L] eq 1 then
				newij := [ij[1]+1,ij[2],
					ij[3],ij[4],ij[5]+ij[3],ij[6]+ij[4]];
				L := L[1..#L-1];
				n := L[#L];
				n -:= 1;
				L[#L] := n;
				n := R[#R];
				n -:= 1;
				R[#R] := n;
			elif (L[#L] eq 1 and R[#R] eq 1 and #R eq 3) or
							R[#R] eq 1 then
				newij := [ij[1],ij[2]+1,
					ij[3]+ij[5],ij[4]+ij[6],ij[5],ij[6]];
				R := R[1..#R-1];
				n := L[#L];
				n -:= 1;
				L[#L] := n;
				n := R[#R];
				n -:= 1;
				R[#R] := n;
			end if;
			Append(~result,newij);
		end while;
		return Reversion(result);
	else
		result := [ [#L-1,#R-1,1,0,0,1] ];
		while #L gt 2 or #R gt 2 do
			ij := result[#result];
			if (L[1] eq 1 and R[1] eq 1 and #L eq 3) or
							L[1] eq 1 then
				newij := [ij[1]-1,ij[2],
					ij[3],ij[4],ij[5]+ij[3],ij[6]+ij[4]];
				L := L[2..#L];
				n := L[1];
				n -:= 1;
				L[1] := n;
				n := R[1];
				n -:= 1;
				R[1] := n;
			elif (L[1] eq 1 and R[1] eq 1 and #R eq 3) or
							R[1] eq 1 then
				newij := [ij[1],ij[2]-1,
					ij[3]+ij[5],ij[4]+ij[6],ij[5],ij[6]];
				R := R[2..#R];
				n := L[1];
				n -:= 1;
				L[1] := n;
				n := R[1];
				n -:= 1;
				R[1] := n;
			end if;
			Append(~result,newij);
		end while;
		return result;
	end if;
end function;

