//////////////////////////////////////////////////////////////////////
//		Diptych varieties. II: Polar varieties
//
// G. Brown, M. Reid.
// August 2012
//
// Magma code for the calculations of the paper
//		Diptych varieties. II: Polar varieties
// available from
// http://www-staff.lboro.ac.uk/~magdb/aflip.html
//
// Contents (section numbers refer to the paper):
//	3.3		[3,1,3]
//	3.4 Part 1	[3,1,3,1]
//	3.4 Part 2	[1,3,1,3]
//	3.4 Observation	[1,3,1,3] with extra annotation
//
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
// Section 3.3
//			[3,1,3]
//////////////////////////////////////////////////////////////////////

RR<x0,x1,x2,x3,y0,y1,y2,A,B,L,M> := PolynomialRing(Rationals(),11);
// The two equations x1*y0 = ... and x0*y1 = ... of the bottom cross are given.
// First compute the ideal of the diptych using usual pentagram method.
// Pentagram x1,x0,y0,y1 to add x3:
M1 := -AntisymmetricMatrix([ x1, B,x0, -L,A*B,y0, -x3,-M*y0,x1,y1 ]);
// Pentagram x1,y0,y1,x3 to add x2:
M2 := -AntisymmetricMatrix([ x1, x3*A*B,y0, -L*M,1,y1, -x2,-B,x1*L,x3 ]);
// Pentagram x3,x1,y0,y1 to add y2:
M3 := -AntisymmetricMatrix([ x3, y1,x1, -L,B,y0, -y2,-L*M,x3*A,y1 ]);
I0 := Ideal(&cat[Pfaffians(m,4) : m in [M1,M2,M3]]);
I0dip := ColonIdeal(I0,y0);
assert IsPrime(I0dip);
// To see the difference - the long equation y0*y2 = ... - compare
MinimalBasis(I0);
MinimalBasis(I0dip);
/* Respectively these are
[
    x3*y1 - y2*B - L^2*M,
    -x1*L + x3*y0 - y1*B,
    -x1*x3 + x2 + B*L*M,
    x1*y2 - x3^2*A - y1*L*M,
    x3*A*L - y0*y2 + y1^2,
    x1*y1 - x3*A*B - y0*L*M,
    -x0*L + x1*y0 - A*B^2,
    -x0*x3 + x1^2 + y0*B*M,
    x0*y1 - x1*A*B - y0^2*M
]
and
[
    x3*y1 - y2*B - L^2*M,
    -x1*L + x3*y0 - y1*B,
    -x1*x3 + x2 + B*L*M,
    x1*y2 - x3^2*A - y1*L*M,
    x3*A*L - y0*y2 + y1^2,
    x1*y1 - x3*A*B - y0*L*M,
    -x0*L + x1*y0 - A*B^2,
    -x0*x3 + x1^2 + y0*B*M,
    x0*y2 - x1*x3*A - y0*y1*M - A*B*L*M,
    x0*y1 - x1*A*B - y0^2*M
]
*/

// And now in a big 6x6 matrix after eliminating x0:
M6 := -AntisymmetricMatrix([ x1, B,x2, -1,L*M,x3, -y0,-x3*A*B,x1*L,y1,
	-y1,-x3^2*A,L^2*M,y2,x3*A*L]);
I6 := Ideal(Pfaffians(M6,4));
assert IsPrime(I6);
assert I6 subset I0dip;

// Reintroducing x0 into a 7x7 matrix:
M7 := -AntisymmetricMatrix([ x1, B,x2, -1,L*M,x3, -y0,-x3*A*B,x1*L,y1,
-y1,-x3^2*A,L^2*M,y2,x3*A*L, x0,A*B^2*M,-x1^2,-M*y0,A*B*x1,x1*x3*A + A*B*L*M]);
I7 := Ideal(Pfaffians(M7,4));
// The following line is useful when amending signs to this correct state:
//	[ f : f in Basis(I7) | not f in I0dip ];
assert ColonIdeal(I7,M) eq I0dip;

// And finally altogether in an extrasymmetric 6x6 matrix including x0:
//	THINK: not yet true extrasymmetric
M8 := -AntisymmetricMatrix([ A*B^2, y0,L, x0,x1,1, x1^2,x2,x3,B*M,
	x1*x3*A+A*B*L*M,x3^2*A,y2,y1*M,L^2*M^2]);
I8 := Ideal(Pfaffians(M8,4));
assert ColonIdeal(I8,L*M) eq I0dip; 


//////////////////////////////////////////////////////////////////////
// Section 3.4, part 1
//			[3,1,3,1]
//////////////////////////////////////////////////////////////////////

RR<x0,x1,x2,x3,x4,y0,y1,y2,A,B,L,M> := PolynomialRing(Rationals(),12);
M1 := -AntisymmetricMatrix([ x4, y1^2,x0, -L,B,y0, -y2,-M,A,y1 ]);
M2 := -AntisymmetricMatrix([ x0, A*B,y0, -M,1,y1, -x2,-y1*B,x0*L,x4 ]);
M3 := -AntisymmetricMatrix([ x0, x2,y1, -B*M,1,x4, -x1,-A*B,x0*L*M,x2 ]);
M4 := -AntisymmetricMatrix([ x2, x4*A*B,x0, -L*M,1,y1, -x3,-B*M,x2,x4 ]);
I := Ideal(&cat[Pfaffians(m,4) : m in [M1,M2,M3,M4]]);
assert IsPrime(I);	// So we already have an irreducible variety
MinimalBasis(I);

/* the minimal basis is:
[
    x0*y1 - y0*M - A*B,
    x0*y2 - x4*A - y1^2*M,
    -y0*y2 + y1^3 + A*L,
    x4*y1 - y2*B - L*M,
    -x0*L + x4*y0 - y1^2*B,
    -x0*x4 + x2 + y1*B*M,
    -x0^2*x4 + x1 + y0*B*M^2 + 2*A*B^2*M,
    -x0*x4^2 + x3 + y2*B^2*M + 2*B*L*M^2
]
*/

//////////////////////////////////////////////////////////////////////
// Section 3.4, part 2
//			[1,3,1,3]
//////////////////////////////////////////////////////////////////////

RR<x0,x1,x2,x3,x4,y0,y1,A,B,L,M> := PolynomialRing(Rationals(),11);
M1 := -AntisymmetricMatrix([ x1, A*B,x0, -L,A,y0, -x3,-M*y0,x1,y1 ]);
M2 := -AntisymmetricMatrix([ x3, y1*B,x1, -L,A,y0, -x4,-L*M,x3,y1 ]);
I := Ideal(&cat[Pfaffians(m,4) : m in [M1,M2]]);
Q,P := PrimaryDecomposition(I);
assert &and[ Q[i] eq P[i] : i in [1..#Q] ]; // primes = primaries, so I reduced
Q[2];	// the residual bit
	// Homogeneous, Dimension 5, Radical, Prime
	// [ x1, x3, y0, y1, A, L ]	--- Note: x2 still in the ring
// Q[1];	// the diptych
assert Q[1] eq ColonIdeal(I,y0); // So it is easy to get Q[1] without PrimDec

M3 := -AntisymmetricMatrix([ x1, x3*A,y0, -L*M,1,y1, -x2,-A*B,x1*L,x3 ]);
Idip := Ideal(Basis(Q[1]) cat Pfaffians(M3,4));

M0 := Module(RR,1);
M1 := quo< M0 | [ [f] : f in Basis(Idip) ] >;
FreeResolution(M1);
//	Chain complex with terms of degree 6 down to -1
//	Dimensions of terms: 0 1 10 25 25 10 1 0

// Observation: [1,3,1,3] with extra annotation C
RR<x0,x1,x2,x3,x4,y0,y1,A,B,C,L,M>
  := PolynomialRing(Rationals(),12);
M1 := -AntisymmetricMatrix([ x1, A*B,x0, -L,A,y0, -x3,-M*y0,x1,y1 ]);
M2 := -AntisymmetricMatrix([ x3, y1*B,x1, -L,A,y0, -x4,-L*M,x3,y1 ]);
I := Ideal(&cat[Pfaffians(m,4) : m in [M1,M2]]);
Q,P := PrimaryDecomposition(I);
assert &and[ Q[i] eq P[i] : i in [1..#Q] ]; // primes = primaries, so I reduced
Q[2];	// the residual bit
	// Homogeneous, Dimension 5, Radical, Prime
	// [ x1, x3, y0, y1, A, L ]	--- Note: x2 still in the ring
// Q[1];	// the diptych
assert Q[1] eq ColonIdeal(I,y0);	// trivial to recover without PrimDec

M3 := -AntisymmetricMatrix([ x1, x3*A,y0, -L*M,C,y1, -x2,-A*B,x1*L,x3 ]);
Idip := Ideal(Basis(Q[1]) cat Pfaffians(M3,4));

time J := ColonIdeal(Idip,y0*L*M);
assert IsPrime(J);

M0 := Module(RR,1);
M1 := quo< M0 | [ [f] : f in Basis(J) ] >;
FreeResolution(M1);
//	Chain complex with terms of degree 7 down to -1
//	Dimensions of terms: 0 1 11 35 50 35 11 1 0

