blowequiv.lib

version="1.0";
category="Miscellaneous";
info="
LIBRARY:   blow-equiv.lib  Fusion of equivariant matrix factorisations
AUTHOR:    Nils Carqueville, Daniel Murfet
KEYWORDS:  matrix factorisation

PROCEDURES:
";

LIB "linalg.lib";
LIB "matrix.lib";
LIB "ring.lib";
LIB "control.lib";
LIB "blownew.lib";

////////////////////////////////////////////////////////////////////
// GROUPS
//
// We assume that the groups involved are finite and abelian, and thus
// may be written as G = Z/(d1) x ... x Z/(dk) for some k \ge 1 where
// d1, ... , dk are integers > 1. So a group is an intvec (d1,...,d_k)
// where k >= 1.

////////////////////////////////////////////////////////////////////
// GROUP REPRESENTATIONS
//
// We specify a G-representation via the equivalence between G-representations and
// graded vector spaces. Recall that a Z/(d) representation splits as a sum of eigenspaces
// for the generator g = 1 of Z/(d) with eigenvalues given by the dth roots of unity.
// So with the primitive root xi chosen, we can specify the structure of a G-representation
// on an m-dimensional space by giving a tuple of non-negative integers (a1,...,am) where
// the ith basis element is acted on by g via multiplication by xi^{ai}. 
//
// To specify a representation of G = Z/(d1) x ... x Z/(dk) on an m-dimensional space
// (note that all the actions can be simultaneously diagonalised) is to specify a list
// (A1, ... , Am) where each Ai = (ai1, ..., aik) is an intvec, specifying the ith
// basis element as a representation of G.

////////////////////////////////////////////////////////////////////
// EQUIVARIANT MATRIX FACTORISATIONS
//
// We expect our ambient ring to be of one of the type
//
//		(0,r),((x(1..),y(1..)),dp
// 
// The x-variables are the "internal" variables. We expect "minpolyblow" to be specified
// see "blow.lib" for instructions.
//
// We assume that G acts linearly on QQ(r)[x1,...,xt] in such a way that the variables
// x(j) are all one-dimensional representations. We may therefore summarise the action
// by specifying the G-grading of this one-dimensional representation, by giving a list
// "varGrading" whose jth entry is an intvec giving the grading of x(j).

////////////////////////////////////////////////////////////////////
// fuseEquivDefects
//
// Returns the reduced tensor product of Y with X over G where Y is
// a G-equivariant matrix factorisation of V(y1) - W(x) and X is a
// G-equivariant matrix factorisation of W(x) - U(y2). Here the group
// G acts only on the x variables, but we assume there may be other
// (abelian) groups acting on V, U and we keep track of these gradings
// as well during the fusion.
//
// The G-equivariant structure on both Y and X is specified by two
// lists of intvecs Yequiv, Xequiv, as specified above where these
// give in order X = X0 + X1 a basis and the G-action on that basis.
//
// Inputs:
//    list roots_of_unity: if G is an intvec (d1, ..., dk) then this list
//    should contain polynomials (xi1, ..., xik) where xi1 is a chosen
//    dith primitive root of unity.
//
// Optional inputs:
//    1. left group F
//    2. F-grading for all y variables
//    3. additional grading Yequiv_extra
//    4. right group H
//    5. H-grading for all y variables
//    6. additional grading Xequiv_extra
//
// Recall that X: (W, H) --> (V, G) and Y: (V, G) --> (U, F). You must specify either
// no additional inputs, 1-3 or 1-6.
//
// Returns: 
//    - if no optional inputs are specified, just returns the matrix of the fusion
//    - if additional inputs are specified, returns a list consisting of
//       - the differential
//       - the left equivariant structure
//       - the right equivariant structure

proc fuseEquivDefects(matrix Y, matrix X, intvec G, list varGrading, list roots_of_unity, poly W, list Yequiv, list Xequiv, list #)
{
	int prl = printlevel;
	
	// Deal with additional optional gradings
	if( size(#) > 0 )
	{
		intvec F = #[1];
		list varGradingF = #[2];
		list Yequiv_extra = #[3];
	}
	
	if( size(#) > 3 )
	{
		intvec H = #[4];
		list varGradingH = #[5];
		list Xequiv_extra = #[6];
	}
	
	if( !isEquivMFValid(Y, Yequiv, G, varGrading) )
	{
		print("[fuseEquivDefects] Passed invalid equivariant structure on Y.");
		return();
	}
	
	if( !isEquivMFValid(X, Xequiv, G, varGrading) )
	{
		print("[fuseEquivDefects] Passed invalid equivariant structure on X.");
		return();
	}
	
	// Note that the theorem says that splitting the idempotent e on (Y|X)^G
	// computes Y *_{A_G} X \theta_1 ... \theta_m where X \theta_1 ... \theta_m
	// denotes the G-equivariant matrix factorisation obtained by tensoring X
	// with the finite-dimensioanl Z2-graded G-representation k \theta_1 ... \theta_m
	// which is a subrepresentation of \wedge F_\theta where F_\theta
	// is the direct sum k \theta_1 \oplus ... \oplus k \theta_m where
	// g * \theta_i = \alpha_{g,i}^{-1} \theta_i where g * x(i) = \alpha_{g_i} x(i)
	// and these \alpha_{g_i} are what varGrading records.
	
	// Hence, in order for this routine to return Y *_{A_G} X we should first
	// take X and perform the opposite twist, i.e. replace X by the G-equivariant
	// MF X \psi_1 ... \psi_m where g * \psi_i = \alpha_{g,i} \psi_i and \psi_i
	// is odd in the Z2-grading.
	
	int numX = numXVars();
	
	// Shift the G-grading
	int i;
	intvec tensor_grade;
	for(i=1; i<=numX; i++)
	{
		tensor_grade = tensor_grade + varGrading[i];
	}
	
	for(i=1; i<=size(Xequiv); i++)
	{
		Xequiv[i] = Xequiv[i] + tensor_grade;
	}
	
	// Shift the Z2-grading
	X = mfSuspend(X, numX);
	Xequiv = mfSuspendEquiv(Xequiv, numX);
	
	if( size(#) > 3 )
	{
		Xequiv_extra = mfSuspendEquiv(Xequiv_extra, numX);
	}
	
	// Construct the jacobi ideal
	ideal I;
	
	for(i=1; i<=numX; i++)
	{
		I = I + ideal( diff(W, x(i)));
	}
	
	// Compute the ordinary tensor product:
	dbprint(printlevel, "[fuseEquivDefects] Inflating tensor product");
	matrix D = MFtensor(Y,X);
	matrix Dinf = inflate(D, I);

	list XYequiv = MFtensorEquiv(Yequiv, Xequiv);
	list Dequiv = inflateEquiv(XYequiv, I, varGrading);
	
	if( size(#) > 0 )
	{
		list Yequiv_zeros;
		list varGrading_zeros;
		list Xequiv_zeros;
		for(i=1; i<= size(Xequiv); i++){ Xequiv_zeros[i] = 0 * Xequiv[i]; }
		for(i=1; i<= size(Yequiv); i++){ Yequiv_zeros[i] = 0 * Yequiv[i]; }
		for(i=1; i<= size(varGrading); i++){ varGrading_zeros[i] = 0 * varGrading[i]; }
		
		list left_equiv = MFtensorEquiv(Yequiv_extra, Xequiv_zeros);
		left_equiv = inflateEquiv(left_equiv, I, varGrading_zeros);
	}
	
	if( size(#) > 3 )
	{
		list right_equiv = MFtensorEquiv(Yequiv_zeros, Xequiv_extra);
		right_equiv = inflateEquiv(right_equiv, I, varGrading_zeros);
	}
	
	// Compute the symmetrised homotopies acting on
	// the matrix factorisation X, as in (20.1) of
	// (orbfus)
	list lambdas;
	
	int i;
	for(i=1; i<=numX; i++)
	{
		list Gelts = groupElements(G);
		matrix A[nrows(X)][ncols(X)];
		
		// \Lambda_i = 1/|G| \Sum_{g \in G} \alpha^{-1}_{g, i} A_g ( g^{-1} * \partial_{x_i}(d_X)) A_g^{-1}
		// where g * x(i) = \alpha_{g,i} x(i)
		
		int j;
		for(j=1; j <= size(Gelts); j++)
		{
			intvec g = Gelts[j];
			
			// s is \alpha_{g^{-1},i}
			poly s = 1;
			int a;
			for(a=1; a <= size(G); a++)
			{
				s = s * roots_of_unity[a]^((G[a] - g[a]) * varGrading[i][a]);
			}
				
			matrix Ag = gmatrix(G, G - g, Xequiv, roots_of_unity);
			matrix Ag_inv = gmatrix(G, g, Xequiv, roots_of_unity);
			matrix diff_act = actmatrix(G, G - g, varGrading, diff(X, x(i)), roots_of_unity);
			
			matrix B = s * Ag * diff_act * Ag_inv;
			A = A + B;
		}
		
		A = 1/poly(size(Gelts)) * A;
		
		lambdas[i] = ZZtensor( unitmat(ncols(Y)), A );
		
		// Check this is a homotopy
		if( A * X + X * A != diff(W, x(i)) * unitmat(nrows(X)) )
		{
			print("[fuseEquivDefects] Failed to construct homotopy.");
			return();
		}
	}
	
	// Construct the idempotent E
	dbprint(printlevel, "[fuseEquivDefects] Constructing idempotent");
	matrix AtProduct = unitmat(nrows(Dinf));
	matrix LambdaProduct = unitmat(nrows(Dinf));
	
	int i;
	matrix z[nrows(Dinf)][ncols(Dinf)];
	for(i=1; i<=numX; i++)
	{
		matrix at = atiyahClass(I, D, i);
		matrix la = inflate( lambdas[i], I );
		
		if( la * Dinf + Dinf * la != z )
		{
			print("[fuseEquivDefects] Lambda " + string(i) + " is not an odd closed map.");
			return();
		}
		
		if( at * Dinf + Dinf * at != z )
		{
			print(at);
			print("[fuseEquivDefects] Atiyah class " + string(i) + " is not an odd closed map.");
			return();
		}
		
		AtProduct = AtProduct * at;
		LambdaProduct = LambdaProduct * la;
	}
	
	int sign = (numX * (numX+1)) div 2;
	matrix E = (-1)^(sign) * LambdaProduct * AtProduct;
	
	if( E * Dinf - Dinf * E != z )
	{
		print("[fuseEquivDefects] E is not a morphism, exiting.");
		return();	
	}
	
	// Compute the projector onto G-degree 0 of the inflation (i.e. G-invariants)
	dbprint(printlevel, "[fuseEquivDefects] Projecting onto G-degree zero");
	list PiG_list;
	
	for(i = 1; i <= size(Dequiv); i++)
	{
		if( modintvec(Dequiv[i],G) == 0 )
		{
			PiG_list = PiG_list + list(i);
		}
	}
	
	if( size(PiG_list) == 0 )
	{
		print("[fuseEquivDefects] G-invariant subspace is zero, so the fusion is the zero matrix factorisation.");
		return(0);
	}
	
	matrix PiG[size(PiG_list)][ncols(Dinf)];
	matrix IiG[ncols(Dinf)][size(PiG_list)];
	
	for(i = 1; i <= size(PiG_list); i++)
	{
		PiG[i,PiG_list[i]] = 1;
		IiG[PiG_list[i],i] = 1;
	}
	
	// The restriction of the differential Dinf to the G-invariant part
	matrix Dinf_G = PiG * Dinf * IiG;
	
	matrix z[nrows(Dinf_G)][ncols(Dinf_G)];
	if( Dinf_G * Dinf_G != (D*D)[1,1] * unitmat(nrows(Dinf_G)) )
	{
		print("[fuseEquivDefects] The G-invariant part is not a matrix factorisation.");
		return();	
	}
	
	// The projection and inclusions should be morphisms of matrix factorisations
	if( Dinf_G * PiG != PiG * Dinf )
	{
		print("[fuseEquivDefects] The projection is not a morphism.");
		return();
	}
	
	if( IiG * Dinf_G != Dinf * IiG )
	{
		print("[fuseEquivDefects] The injection is not a morphism.");
		return();
	}
	
	// Induce gradings for the remaining groups on the G-invariant part
	if( size(#) > 0 )
	{
		left_equiv = computeEquivFromInjectionYVar(IiG, F, left_equiv, varGradingF);
	}
	
	if( size(#) > 3 )
	{
		right_equiv = computeEquivFromInjectionYVar(IiG, H, right_equiv, varGradingH); 
	}
	
	// Reduce this inflated differential:
	dbprint(printlevel, "[fuseEquivDefects] Reducing inflated differential of size " + string(ncols(Dinf_G)));
	list l = mfReduce(Dinf_G);
	matrix RD = l[1];
	matrix sF = l[2]; // A homotopy equivalence sF: Dinf_G -> RD
	matrix sG = l[3]; // A homotopy equivalence sG: RD -> Dinf_G
	
	// Induce gradings
	if( size(#) > 0 )
	{
		left_equiv = computeEquivFromInjectionYVar(sG, F, left_equiv, varGradingF);
	}
	
	if( size(#) > 3 )
	{
		right_equiv = computeEquivFromInjectionYVar(sG, H, right_equiv, varGradingH); 
	}
	
	// The induced homotopy idempotent on the G-invariant part
	matrix E_G = PiG * E * IiG;
	
	if( E_G * Dinf_G - Dinf_G * E_G != z )
	{
		print("[fuseEquivDefects] The induced idempotent on the G-invariant part is not a morphism");
		return();
	}
	
	// The idempotent on the reduction of the G-invariant part
	matrix ep = sF * E_G * sG;
	
	// Sanity check: ep must be an endomorphism of RD:
	int h = ncols(ep * RD);
	matrix zero[h][h];
	
	if( ep * RD - RD * ep != zero )
	{
		print("[fuseEquivDefects] The reduced idempotent is not a morphism, exiting.");
		return();
	}
	
	matrix ep_strict = ep;
	
	if( ep * ep != ep )
	{
		// We need to strictify
		dbprint(printlevel, "[fuseEquivDefects] Strictifying idempotent");
		ep_strict = mfStrictifyIdempotent(RD, ep);
	}
	
	if( ep_strict == zero )
	{
		print("[fuseEquivDefects] The fusion is the zero matrix factorisation.");
		return(ep_strict);	
	}
	
	// Split the idempotent
	dbprint(printlevel, "[fuseEquivDefects] Splitting idempotent");
	list l = mfSplitIdempotent(RD, ep_strict);
	
	if( typeof(l[1]) == "string" )
	{
		print("[fuseEquivDefects] Obtained a Z2-graded complex with different ranks in even and odd degrees. Currently not handled, sorry.");
		return();
	}
	
	matrix A = l[1];
	matrix sG = l[3]; // The split mono
	
	if( size(#) > 3 )
	{
		left_equiv = computeEquivFromInjectionYVar(sG, F, left_equiv, varGradingF);
		right_equiv = computeEquivFromInjectionYVar(sG, H, right_equiv, varGradingH); 
		list l = A, left_equiv, right_equiv;
		return(l);
	}
	
	if( size(#) > 0 )
	{
		left_equiv = computeEquivFromInjectionYVar(sG, F, left_equiv, varGradingF);
		
		list l = A, left_equiv;
		return(l);
	}
	
	return(A);
}

proc modintvec(intvec g, intvec h)
{
	intvec m;
	int i;
	for(i=1; i <= size(g); i++)
	{
		m[i] = g[i] mod h[i];
	}
	
	return(m);
}

proc mfSuspendEquiv(list g, list #)
{
	int power = 1;
	
	if( size(#) != 0 ){ power = #[1]; }
	
	if( power mod 2 == 0 )
	{
		return(g);
	}
	else
	{
		int n = size(g) div 2;
		list r = g[(n+1)..2*n], g[1..n];
	
		return(r);
	}
}

////////////////////////////////////////////////////////////////////
// actpoly
//
// Implements the action of G on the polynomial ring
//
// Inputs:
//     G = (d1, ..., dk)
//     g = (a1, ..., ak) gives a group element
//     varGrading is a list specifying the representation of G on x variables
//     polynomial f

proc actpoly(intvec G, intvec g, list varGrading, poly f, list roots_of_unity)
{
	int numX = numXVars();
	
	int i, a;
	poly r = f;
	
	for(i=1; i <= numX; i++)
	{
		// s is \alpha_{g,i}
		poly s = 1;
		for(a=1; a <= size(G); a++)
		{
			s = s * roots_of_unity[a]^(g[a] * varGrading[i][a]);
		}
		
		r = subst(f, x(i), s * x(i));
	}
	
	return(r);
}

proc actmatrix(intvec G, intvec g, list varGrading, matrix F, list roots_of_unity)
{
	matrix R[nrows(F)][ncols(F)];
	
	int i,j;
	for(i=1; i<=nrows(F); i++)
	{
		for(j=1; j<=ncols(F); j++)
		{
			R[i,j] = actpoly(G, g, varGrading, F[i,j], roots_of_unity);
		}
	}
	
	return(R);
}


////////////////////////////////////////////////////////////////////
// gmatrix
//
// Inputs:
//    intvec G = (d1, ..., dk)
//    intvec g = (a1, ..., ak) with 0 \le ai > di
//    Xequiv = grading matrix for a G-representation
//    roots_of_unity = (xi1, ..., xik)
//
// Returns:
//    matrix of g acting on X

proc gmatrix(intvec G, intvec g, list Xequiv, list roots_of_unity)
{
	int ra = size(Xequiv);
	matrix Ag[ra][ra];
	
	int i, a;
	for(i=1; i<=ra; i++)
	{
		intvec y = Xequiv[i]; // y = (y1,...,yk)
		
		// G = Z/(d1) x ... x Z/(dk) and the generator
		// of Z/(dj) acts on the ith generator of X by
		// multiplication with xi(j)^{yj}. Hence the element
		// g acts on this generator by multiplication with
		//
		//   xi(1)^{a1y1} ... xi(k)^{akyk}
		 
		Ag[i,i] = 1;
		
		for(a=1; a<=size(G); a++)
		{
			Ag[i,i] = Ag[i,i] * roots_of_unity[a]^(g[a] * y[a]);
		}
	}
	
	return(Ag);
}


////////////////////////////////////////////////////////////////////
// groupElements
//
// Inputs:
//    intvec G = (d1, ..., dk)
//
// Returns the list of all intvecs (a1, ..., ak) where 0 \le ai < di
// for all 1 \le i \le k.

proc groupElements(intvec G)
{
	if( size(G) == 1 )
	{
		list elts;
		int i;
		for(i=0; i < G[1]; i++)
		{
			intvec t = i;
			elts = elts + list(t);
		}
		
		return(elts);
	}
	
	intvec subG = G[2..size(G)];
	list sublist = groupElements(subG);
	
	list elts;
	int i,j;
	for(i=0; i < G[1]; i++)
	{
		for(j=1; j<=size(sublist); j++)
		{
			intvec t = i, sublist[j];	
			elts = elts + list(t);
		}
	}
	
	return(elts);
}

////////////////////////////////////////////////////////////////////
// MFtensorEquiv
//
// Returns the equivariant structure on the tensor product of two
// equivariant matrix factorisations

proc MFtensorEquiv(list g, list h)
{
	// With the notation of MFtensor, X x Y is, as a Z/2-graded free module
	// (X x Y)0 = (X0 x Y0) + (X1 x Y1)
	// (X x Y)1 = (X0 x Y1) + (X1 x Y0)
	// with the basis ordered in the usual way. If X,Y are free modules with G-representations
	// with the gradings described by lists g, h resp. then the intvec describing
	// the G-representation on the free modules in X x Y is
   
	int sizeX = size(g) div 2;
	int sizeY = size(h) div 2;
	int rankProduct = sizeX * sizeY;
	
	list r;
	
	int i, j;
	for(i=1;i<=sizeX;i++)
	{
		for(j=1;j<=sizeY;j++)
		{
			r[(i-1)*sizeY + j] = g[i] + h[j]; // Gradings on X0 x Y0
			r[rankProduct + (i-1)*sizeY + j] = g[sizeX + i] + h[sizeY + j]; // Grading on X1 x Y1
			r[2*rankProduct + (i-1)*sizeY + j] = g[i] + h[sizeY + j]; // X0 x Y1
			r[3*rankProduct + (i-1)*sizeY + j] = g[sizeX + i] + h[j]; // X1 x Y0
		}
	}
	
	return(r);
}

//////////////////////////////////////////////////////////////////////////////////
// "inflateEquiv" returns the G-grading vector for an inflated MF, given the G-grading
// vector of the original MF and the ideal to inflate with. We use the basis for 
// R/J given by reduce(kbase(std(J)),std(J)).
//
// The underlying graded free module of the MF is R{g[1]} + R{g[2]} + ... and
// the underlying graded free module of the inflation is a direct sum of modules
// of the form R{g[i]} x R/J which has the usual grading.
//////////////////////////////////////////////////////////////////////////////////

proc inflateEquiv(list g, ideal J, list varGrading)
{
	def RRR = basering;
	def nR = ringWithoutYVars();
	setring nR;
		
	ideal J = imap(RRR,J);
	ideal Jstd = std(J);
	ideal basis = kbase(Jstd);
	int jdim = size(basis);
	module BB = reduce(basis,Jstd);
	int numX = numXVars();
	
	list r;   
	
	int i,j,k;
	for(i=1;i <= size(g); i++)
	{
		for(j=1;j<=jdim;j++)
		{
			// The basis element BB[j] is a monomial in the x variables
			// and each x variable x(i) is a representation with G-grading
			// varGrading[i]. Hence the G-grading of x(1)^a1 ... x(t)^at
			// is the intvec a1 * varGrading[1] + ... + at * varGrading[t]
			
			intvec gengrade;
			for(k=1; k<=numX; k++)
			{
				gengrade = gengrade + leadexp(BB[j])[k] * varGrading[k];
			}
			
			//print("combination of " + string(g[i]) + " with " + string(BB[j]));
			r[(i-1)*jdim + j] = g[i] + gengrade;
		}
	}
	
	setring RRR;
	
	return(r);
}

////////////////////////////////////////////////////////////////////
// computeEquivFromInjection
//
// Given a matrix F interpreted as a morphism of free modules, and a
// list g assigning equivariant structure to the basis elements of the target module,
// we derive the induced grading intvec of the source module making the map F
// into a homogeneous map of degree zero. Note that these gradings should be treated
// modulo the order of the group.
//
// Inputs:
//     g is a list of size nrows(F), with g[j] giving the G-grading of
//     the jth generator, meaning that g[j] is an intvec of size |G|.

proc computeEquivFromInjection(matrix F, intvec G, list g, list varGrading)
{   
	// Check that the sizes of g and F are compatible:
	if( size(g) != nrows(F) )
	{
		print("[computeEquivFromInjection] Sizes of map and grading vector not compatible, exiting.");
		return();
	}
	
	if( F == matrix(0,nrows(F),ncols(F)) )
	{
		print("[computeEquivFromInjection] The given map is zero, exiting.");
		return();
	}
	
	int i,j,k;
	list r;
	int ncolsF = ncols(F);
	int sizeg = size(g);
	
	list varNames = ringlist(basering)[2];
	
	int a;
	poly product_x_vars = 1;
	
	for(a = 1; a <= size(varNames); a++)
	{
		if( varNames[a][1] == "x" )
		{
			product_x_vars = product_x_vars * var(a);
		}
	}
	
	int b;
	
	for(i=1;i<=ncolsF;i++)
	{
		list rowDegrees;
		
		for(j=1;j<=sizeg;j++)
		{
			if( F[j,i] != 0 )
			{
				//print("nonzero entry F[" + string(j) + "," + string(i) + "] = " + string(F[j,i]));
				
				// First row of M are monomials in x variables
				// second row are the coefficients (in y variables)
				matrix M = coef( F[j,i], product_x_vars );
				
				// Each monomial that appears should have a consistent G-grade
				for(b=2; b <= ncols(M); b++)
				{
					if( modintvec( gradeMonom(M[1,b], varGrading, "x") - gradeMonom(M[1,b-1], varGrading, "x"), G ) != 0 )
					{
						print("[computeEquivFromInjection] monomials in entry F[" + string(j) + "," + string(i) + "] = " + string(F[j,i]) + " do not have consistent G-degree.");
						return();
					}
				}
				
				intvec d = modintvec( g[j] + gradeMonom(M[1,1], varGrading, "x"), G );
				
				// Here d is the predicted degree of the ith generator of the
				// domain free module, according to the entry F[j,i]
				rowDegrees = rowDegrees + list(d);
			}
		}
		
		// Check that each nonzero entry in this column agrees on the degree
		for(k=2; k <= size(rowDegrees); k++)
		{
			if( modintvec( rowDegrees[k] - rowDegrees[k-1], G ) != 0 )
			{
				print("[computeEquivFromInjection] Encountered multiple nonzero degrees in column i = " + string(i) + " exiting.");
				print("[computeEquivFromInjection] Competing degrees are " + string(rowDegrees[k]) + " and " + string(rowDegrees[k-1]));
				return();
			}
		}
		
		r[i] = rowDegrees[1];
	}
	
	return(r);
}

proc gradeMonom(poly m, list varGrading, string v)
{
	list varNames = ringlist(basering)[2];
	
	intvec gengrade;
	int i, j;
	for(i=1; i<=size(varNames); i++)
	{
		if( varNames[i][1] == v )
		{
			j = j + 1;
			
			gengrade = gengrade + leadexp(m)[i] * varGrading[j];
		}
	}

	return(gengrade);
}

////////////////////////////////////////////////////////////////////
// computeEquivFromInjectionYVar
//
// Same as above, but now we are given a group acting on the y variables
// and gradings of those variables.

proc computeEquivFromInjectionYVar(matrix F, intvec G, list g, list varGrading)
{   
	// Check that the sizes of g and F are compatible:
	if( size(g) != nrows(F) )
	{
		print("[computeEquivFromInjectionYVar] Size of grading vector (" + string(size(g)) + ") and number of rows in the map (" + string(nrows(F)) + ") not equal, exiting.");
		return();
	}
	
	if( F == matrix(0,nrows(F),ncols(F)) )
	{
		print("[computeEquivFromInjectionYVar] The given map is zero, exiting.");
		return();
	}
	
	int i,j,k;
	list r;
	int ncolsF = ncols(F);
	int sizeg = size(g);
	
	list varNames = ringlist(basering)[2];
	
	int a;
	poly product_y_vars = 1;
	
	for(a = 1; a <= size(varNames); a++)
	{
		if( varNames[a][1] == "y" )
		{
			product_y_vars = product_y_vars * var(a);
		}
	}
	
	int b;
	
	for(i=1;i<=ncolsF;i++)
	{
		list rowDegrees;
		
		for(j=1;j<=sizeg;j++)
		{
			if( F[j,i] != 0 )
			{
				//print("nonzero entry F[" + string(j) + "," + string(i) + "] = " + string(F[j,i]));
				
				// First row of M are monomials in y variables
				// second row are the coefficients (in x variables)
				matrix M = coef( F[j,i], product_y_vars );
				
				// Each monomial that appears should have a consistent G-grade
				for(b=2; b <= ncols(M); b++)
				{
					if( modintvec( gradeMonom(M[1,b], varGrading, "y") - gradeMonom(M[1,b-1], varGrading, "y"), G ) != 0 )
					{
						print("[computeEquivFromInjectionYVar] monomials in entry F[" + string(j) + "," + string(i) + "] = " + string(F[j,i]) + " do not have consistent G-degree.");
						return();
					}
				}
				
				intvec d = modintvec( g[j] + gradeMonom(M[1,1], varGrading, "y"), G );
				
				// Here d is the predicted degree of the ith generator of the
				// domain free module, according to the entry F[j,i]
				rowDegrees = rowDegrees + list(d);
			}
		}
		
		// Check that each nonzero entry in this column agrees on the degree
		for(k=2; k <= size(rowDegrees); k++)
		{
			if( modintvec( rowDegrees[k] - rowDegrees[k-1], G ) != 0 )
			{
				print("[computeEquivFromInjectionYVar] Encountered multiple nonzero degrees in column i = " + string(i) + " exiting.");
				print("[computeEquivFromInjectionYVar] Competing degrees are " + string(rowDegrees[k]) + " and " + string(rowDegrees[k-1]));
				return();
			}
		}
		
		r[i] = rowDegrees[1];
	}
	
	return(r);
}

////////////////////////////////////////////////////////////////////
// defectBrunnerRoggenkamp
//
// Returns the Brunner Roggenkamp defect P^{(m,n)} from their paper
// "Defects and bulk perturbations of boundary Landau-Ginzburg orbifolds"
// Section 5.1 (5.5), (5.6) as a Q-graded, Z/d'Z x Z/dZ equivariant MF
// of y^{d'} - x^d. Here we are given an integer 0 \le m < d and
// a vector of d' integers n non-negative with \sum_i n_i = d'.
//
// Assumes variables x(1), y(1) are in the ambient ring and that
// we are passed a primitive d'th root of unity.

proc defectBrunnerRoggenkamp(int m, intvec n, poly root_of_unity)
{   
	int dprime = size(n);
	int d = sum(n);
	
	int a,b;
	matrix Xi[dprime][dprime];

	for(a = 1; a <= dprime; a++)
	{
		for(b = 1; b <= dprime; b++)
		{
			if( (a - b - 1) mod dprime == 0 )
			{
				Xi[a,b] = x(1)^n[a];
			}
		}
	}
	
	if( power(Xi,dprime) != x(1)^d * unitmat(dprime) )
	{
	   print("[defectBrunnerRoggenkamp] Construction of Xi is broken.");
	   return();
	}
	
	matrix p1[dprime][dprime] = y(1) * unitmat(dprime) - Xi;
	
	matrix p0[dprime][dprime] = unitmat(dprime);
	
	int i;
	for(i = 1; i <=dprime-1; i++)
	{
		p0 = p0 * (y(1) * unitmat(dprime) - root_of_unity^i * Xi);
	}
	
	matrix z[dprime][dprime];
	matrix Y = blockmat(z,p1,p0,z);
	
	if( Y * Y != (y(1)^dprime - x(1)^d) * unitmat(nrows(Y)) )
	{
		print("[defectBrunnerRoggenkamp] Construction of P is broken.");
		return();
	}
	
	list gradings_left0;
	list gradings_right0;
	list gradings_left1;
	list gradings_right1;
	
	int ncount = -m;
	
	for(i = 0; i <= dprime - 1; i++)
	{
		if( i > 0 ){ ncount = ncount - n[i+1]; }
		gradings_left0 = gradings_left0 + list(intvec(i));
		gradings_right0 = gradings_right0 + list(intvec(ncount));
		gradings_left1 = gradings_left1 + list(intvec(i+1));
		gradings_right1 = gradings_right1 + list(intvec(ncount));
	}
	
	list gradings_left = gradings_left0 + gradings_left1;
	list gradings_right = gradings_right0 + gradings_right1;
	
	list l = Y, gradings_left, gradings_right;
	
	return(l);
}

////////////////////////////////////////////////////////////////////
// boundaryBrunnerRoggenkamp
//
// Returns the Brunner Roggenkamp boundary condition Q^{(M,N)} from their paper
// "Defects and bulk perturbations of boundary Landau-Ginzburg orbifolds"
// Section 5.3 (5.22) as a Q-graded, Z/dZ equivariant MF
// of x^d. Here we are given integers 0 \le M, N < d.
//
// Assumes variable x(1) is in the ambient ring.

proc boundaryBrunnerRoggenkamp(int M, int N, int d)
{   
	matrix X[2][2] = 0, x(1)^N, x(1)^(d-N), 0;
	list Xequiv;
	
	Xequiv[1] = intvec(M); // G-grading of even generator
	Xequiv[2] = intvec(M+N); // G-grading of odd generator
	
	list l = X, Xequiv;
	
	return(l);
}

////////////////////////////////////////////////////////////////////
// isEquivMFValid
//
// Returns True (1) or False (0) depending on whether or not the
// given data constitutes a G-equivariant matrix factorisation of
// the given potential.
//
// Inputs:
//   

proc isEquivMFValid(matrix D, list Dequiv, intvec G, list varGrading)
{
	list Dequiv_inf = computeEquivFromInjection(D, G, Dequiv, varGrading);
	
	int i;
	for(i = 1; i <= size(Dequiv); i++)
	{
		if( modintvec(Dequiv_inf[i] - Dequiv[i], G) != 0 )
		{
			print("[isEquivMFValid] G-grading not valid at index i = " + string(i) + ".");
			return(0);
		}
	}
	
	return(1);
}