moduli.lib

version="1.0";
category="Miscellaneous";
info="
LIBRARY:   moduli.lib  Moduli of matrix factorisations
AUTHOR:    Daniel Murfet
KEYWORDS:  matrix factorisation

PROCEDURES:
";

// NOTE: We include a version of matrix.lib which suppresses some
// unnecessary output from the procedure "rowred"
LIB "linalg.lib";
LIB "matrix.lib";
LIB "ring.lib";

////////////////////////////////////////////////////////////////////
// USAGE GUIDE
//
// Matrix factorisations are encoded as in blow.lib and other libraries,
// that is, as odd supermatrices with polynomial entries. To create such
// a block matrix from a pair of matrices A, B such that A * B = B * A = W . I,
// use mfblock(A,B), see blow.lib.
//
// When forming the moduli spaces, moduli variables for MFs look like
// f(z)(i)(j)(a) or f(i)(j)(a) where z = 1, 2 means respectively the upper
// right or lower left block, i,j are row and column indices and a is an
// index into a list of monomials. To see what this list of monomials is, run
// displayIndices with the input k you pass the other functions.
//
// Moduli variables for homotopies are h(t)(z)(i)(j)(a) where t = 1,...,r
// is the index into the set of variables of the ring.
//
// Graded MFs: our conventions are as explained in blow.lib (see 
// isGradingValid for the gory details).
//
// TODO: We do not handle algebraic extensions

////////////////////////////////////////////////////////////////////
// moduliOfPotential_H
//
// H for homotopies
//
// Given a polynomial ring R, a polynomial W in R, a rank r, a
// a degree cutoff k, and an integer N_x for each variable x, returns
// a quotient ring which parametrises matrix factorisations of
// W of rank r, with only monomials of total degree <= k, and
// for each variable x a homotopy with only monomials of total 
// degree <=k for the action of x^(N_x).
//
// We assume R = basering. The return value is a ring nR and in this ring
// is already declared an ideal moduli_eqns which is the quotient, 
// i.e. nR/moduli_eqns is the parametrising space.

proc moduliOfPotential_H(poly W, int r, int k, intvec N)
{
	if( size(N) != nvars(basering) )
	{
		print("[moduliOfPotential_H] N has wrong size, exiting.");
		return();
	}
	
	// Get the ring nR with the f(z)(i)(j)(a) variables and
	// the moduli_eqns. To this we are going to add h(t)(z)(i)(j)(a) variables
	// for the homotopies, and then impose additional equations in a final
	// ring hR that gets returned
	
	def nR = moduliOfPotential(W,r,k); 

	// Monomials of total weight <= k in the ring variables, as intvecs
	int numVars = nvars(basering);
	list monoms = monomialdict(k, numVars);
 	int Q = size(monoms);
 	
 	list var_exponents;
 	int t;
 	for( t=1; t<=numVars; t++ )
 	{
 		var_exponents = var_exponents + list(N[t] * leadexp(var(t)));
 	}
 	
 	// The variables of our new ring are
    // h(t)(z)(i)(j)(a) from f(1..numVars)(1,2)(1..r)(1..r)(1..Q)
    // where z = 1 means top right block A, z = 2 means bottom left B,
    // i and j are the row and column resp., and a is the monomial index
    list newVarNames;
    int z,i,j,a,t;
    for( t=1; t<=numVars; t++ )
    {
    	for( z=1; z<=2; z++ )
    	{
    		for( i=1; i<=r; i++ )
	    	{
    			for( j=1; j<=r; j++ )
	    		{
    				for( a=1; a<=Q; a++ )
    				{
    					// WARNING: If you change the ordering here you should also change it in e.g.
	    				// universalMF and other places where the order is used
    					string nv = "h(" + string(t) + ")(" + string(z) + ")(" + string(i) + ")(" + string(j) + ")(" + string(a) + ")";
    					newVarNames = newVarNames + list(nv);
	    			}
    			}
    		}
	    }
 	}
 	
 	// Now create a temporary ring with just the h variables
    def RRR = basering;
    
	list tempR_list;
	tempR_list[1] = ringlist(RRR)[1]; // TOOD: minpoly
	tempR_list[2] = newVarNames;
	
	intvec kk = (1..size(newVarNames));
    for(i=1; i<=size(newVarNames); i++)
    {
        kk[i] = 1;
    }
	tempR_list[3] = list( list( "dp", kk ) );
	tempR_list[4] = ideal(0);
	
	def tempR = ring(tempR_list);
 	
	def hR = nR + tempR;
	
	setring(hR);
	
	// imap works to import stuff from RRR and nR to uR
	ideal moduli_eqns = imap(nR, moduli_eqns);
	export(moduli_eqns);
	
	// Let H(t)(z) denote the rxr matrix whose (i,j) entry has coefficient
	// h(t)(z)(i)(j)(a) of the monomial with index a, and let H(t) denote
	// the odd block matrix with blocks H(t)(1), H(t)(2). Then for each t 
	// between 1 and numVars we have an equation
	//
	// H(t) * D + D * H(t) = x(t)^N[t] . I
	//
	// where D is the block matrix [[0,A],[B,0]]. This means that for
	// each t we have two equations
	//
	// H(t)(1) * B + A * H(t)(2) = x(t)^N[t] . I    (e1)
	// H(t)(2) * A + B * H(t)(1) = x(t)^N[t] . I    (e2)
	//
	// This is an equation in the ring containing all the base ring variables,
	// the f's and the h's. So for each (i,j) and for each monomial we get an
	// equation in the f's and h's only. The monomials that can possibly give
	// nontrivial equations are (a) any monomial of total weight <= 2k and 
	// (b) the monomial x(t)^N[t]
	
    // First make a list of such monomials (as intvecs)
    list eqn_monoms = monomialdict(2*k, numVars);
	
    list homotopy_eqn_list;
    int u;

	for( t=1; t<=numVars; t++ )
    {
    	// Add the intvec of the monomial x(t)^N[t]
    	list this_t_eqn_monoms = eqn_monoms + list(var_exponents[t]);
    	
    	for( i=1; i<=r; i++ )
	    {
    		for( j=1; j<=r; j++ )
	    	{	
    			for( u=1; u<=size(this_t_eqn_monoms); u++ )
	    		{
    				poly e1,e2;
	    			intvec m = this_t_eqn_monoms[u];
    			
    				// The contribution of the right hand side. This is
    				// zero unless i = j and the monomial we are looking at
    				// is var_exponents[t]
    				
	    			poly cw = 0;
    				if( i == j && m == var_exponents[t] )
    				{
    					cw = 1;
    				}
    			
    			    e1 = -cw;
    				e2 = -cw;
    			
    				// Now construct the left hand side
    				// H(t)(1) * B + A * H(t)(2) = x(t)^N[t] . I    (e1)
					// H(t)(2) * A + B * H(t)(1) = x(t)^N[t] . I    (e2)
					
	    			int k, ac, ad;
    				for( k=1; k<=r; k++ )
    				{
    					for( ac=1; ac<=Q; ac++ )
    					{
	    					for( ad=1; ad<=Q; ad++ )
    						{
    							if( monoms[ac] + monoms[ad] == m )
    							{
    								e1 = e1 + h(t)(1)(i)(k)(ac) * f(2)(k)(j)(ad);
    								e1 = e1 + f(1)(i)(k)(ac) * h(t)(2)(k)(j)(ad);
    								
    								e2 = e2 + h(t)(2)(i)(k)(ac) * f(1)(k)(j)(ad);
    								e2 = e2 + f(2)(i)(k)(ac) * h(t)(1)(k)(j)(ad);
	    						}
    						}
    					}
    				}
    				
    				homotopy_eqn_list = homotopy_eqn_list + list(e1,e2);
    			}
    		}
    	}
    }
    
    // convert to a row matrix
    matrix mat_eqn[1][size(homotopy_eqn_list)];
    for( i=1; i<=size(homotopy_eqn_list); i++ )
    {
    	mat_eqn[1,i] = homotopy_eqn_list[i];
    }
    
    ideal homotopy_eqns = ideal(mat_eqn);    
    export(homotopy_eqns);
            
	setring RRR;
	return(hR);
}

////////////////////////////////////////////////////////////////////
// universalMF_H
//
// Given the same input data as moduliOfPotential_H returns the
// universal matrix factorisation of that moduli space, together
// with its homotopies. This is defined over the ring R + hR
// where R is the basering and hR is the output of moduliOfPotential_H. 
//
// The return value is a ring uR and in this ring is already
// declared the ideals moduli_eqns,homotopy_eqns from moduliOfPotential_H and
// a matrix factorisation U in supermatrix form (the universal MF).
// Note that this only satisfies U * U = W . I up to moduli_eqns.
// There are homotopies H(1),...,H(r) where r is the number of variables
// in the basering, which satisfies [H(t),U] = x(t)^N[t] up to
// the ideal homotopy_eqns.

proc universalMF_H(poly W, int r, int k, intvec N)
{
	def RRR = basering;
	def hR = moduliOfPotential_H(W,r,k,N);
	def uR = RRR + hR; 
	
	setring(uR);
	
	// imap works to import stuff from RRR and nR to uR
	ideal moduli_eqns = imap(hR, moduli_eqns);
	ideal homotopy_eqns = imap(hR, homotopy_eqns);
	poly W = imap(RRR,W);
	export(W);
	export(moduli_eqns);
	export(homotopy_eqns);
	
	// Now construct the universal matrix U
	matrix A(1)[r][r]; matrix A(2)[r][r];
	
	list varNames = ringlist(RRR)[2];
 	int numVars = size(varNames);
 	
 	int i;
 	list monoms = monomialdict(k, numVars);
	int Q = size(monoms);

	int z,i,j,a;
    for( z=1; z<=2; z++ )
    {
    	for( i=1; i<=r; i++ )
    	{
    		for( j=1; j<=r; j++ )
    		{
    			for( a=1; a<=Q; a++ )
    			{
    				// The RRR variables come first, so if we call
    				// monomial with something of length <= numVars
    				// we get the monomial from RRR
    				A(z)[i,j] = A(z)[i,j] + f(z)(i)(j)(a) * monomial(monoms[a]);
    			}
    		}
    	}
    }
	
	matrix U = mfblocks(A(1),A(2));
	export(U);
	
	// And the homotopies H(1),...,H(numVars)
	matrix H(1..numVars);
	
	int t;
	for( t=1; t<=numVars; t++ )
	{
		A(1) = zeromat(r);
		A(2) = zeromat(r);
		
		for( z=1; z<=2; z++ )
	    {
    		for( i=1; i<=r; i++ )
    		{
    			for( j=1; j<=r; j++ )
	    		{
    				for( a=1; a<=Q; a++ )
    				{
	    				A(z)[i,j] = A(z)[i,j] + h(t)(z)(i)(j)(a) * monomial(monoms[a]);
    				}
    			}
    		}
	    }
	    
	    H(t) = mfblocks(A(1),A(2));
		export(H(t));
	}
		
	setring RRR;
	return(uR);
}

////////////////////////////////////////////////////////////////////
// maxdegForGradingVector
//
// Given Wdeg, r, g as in locusGradedMF, returns the largest
// total degree possible for a monomial in such a graded MF

proc maxdegForGradingVector(int Wdeg, int r, intvec g)
{
	int maxdeg = 0;
	int n = Wdeg - 1;
	
	int i,j,a;
	
	for(i=1; i<=r; i++)
	{
		for(j=1; j<=r; j++)
		{
	    	int P = n + 1 - g[r+i] + g[j];
	    	int Q = n + 1 - g[i] + g[r+j];
	    	
	    	if( P mod 2 == 0 && P > 2 * maxdeg )
	    	{
	    		maxdeg = P div 2;
	    	}
	    	
	    	if( Q mod 2 == 0 && Q > 2 * maxdeg )
	    	{
	    		maxdeg = Q div 2;
	    	}
		}
	}
	
	return(maxdeg);
}

////////////////////////////////////////////////////////////////////
// locusGradedMF
//
// Returns the equations cutting out graded MFs among all MFs. This is
// based on isGradingValid from blow.lib.
//
// IMPORTANT: this procedure should be called with a current basering
// which is the output of e.g. moduliOfPotential, that is, it should
// contain variables f(z)(i)(j)(a) for the same parameters r, k.
//
// We are passed the degree of W (the ordinary degree, not the doubled
// thing discussed in blow.lib), and the intvec encoding the grading
// which should be of size 2 * r

proc locusGradedMF(int numVars, int Wdeg, int r, int k, intvec g)
{
	if( size(g) != 2 * r )
	{
		print("[locusGradedMF] g has wrong size, exiting.");
		return();
	}
	
	if( maxdegForGradingVector(Wdeg, r, g) > k )
	{
		print("[locusGradedMF] Polite warning: with this grading vector the degree cutoff k will truncate the moduli space.");	
	}
	
	int n = Wdeg - 1;
	
	// Monomials of total weight <= k in the ring variables, as intvecs
	list monoms = monomialdict(k, numVars);
 	int Q = size(monoms);
	
	// Pasted from blow.lib:
	// The doubled grading is adopted for the following reason: let R be a polynomial ring
	// with this grading, and W an element of R with deg(W) = n+1. If W = fg for homogeneous f,g
	// then we can write W: R -> R as a composite of two maps g: R -> R{n+1-|g|} and
	// f: R{n+1-|g|} -> R of graded R-modules of degree n+1.

	// Our matrix is in the shape (0 A newline B 0), in the notation of
	// isGradingValid this is (0 a1 newline a0 0).
	
	int i,j,a;
	list grading_eqn_list;
	
	for(i=1; i<=r; i++)
	{
		for(j=1; j<=r; j++)
		{
	    	// To be valid the polynomial in B[i,j] must be homogenous, and it must
		    // define a morphism of graded modules of degree n + 1
       	    // R{a} -> R{b} where a = g[j] and b = g[r + i]
            // that is, we must have
	        // n+1 - g[r+i] + g[j] == 2 * deg(B[i,j])
            // This means that we impose an equation which zeros out any
       	    // f(2)(i)(j)(a) where a indexes a monomial that does NOT satisfy
           	// this identity

			for( a=1; a<=Q; a++ )
    		{
				if( n + 1 - g[r+i] + g[j] != 2 * sum(monoms[a]) )
				{
					grading_eqn_list = grading_eqn_list + list(f(2)(i)(j)(a));
				}			
			
				if( n + 1 - g[i] + g[r+j] != 2 * sum(monoms[a]) )
				{
					grading_eqn_list = grading_eqn_list + list(f(1)(i)(j)(a));
				}
			}
		}
	}
	
	// convert to a row matrix
    matrix mat_eqn[1][size(grading_eqn_list)];
    for( i=1; i<=size(grading_eqn_list); i++ )
    {
    	mat_eqn[1,i] = grading_eqn_list[i];
    }
    
    ideal grading_eqns = ideal(mat_eqn);    
	return(grading_eqns);	
}

////////////////////////////////////////////////////////////////////
// moduliOfPotential
//
// Given a polynomial ring R, a polynomial W in R, a rank r, and
// a degree cutoff k, returns a quotient ring which parametrises
// matrix factorisations of W of rank r, with only monomials of
// total degree <= k.
//
// Ref: p. 3 cutsys45 and p. 6 cutsys46
//
// Details: a matrix factorisation of W is a pair (A,B) of r x r matrices
// where we assume here that r >= 1. Corresponding to each position (i,j) in A
// and monomial in R of total degree <= k we introduce a new variable f(1)(i)(j)(a)
// where a is some index. Similarly we introduce new variables f(2)(i)(j)(a) for B.
// These are the variables of our new ring, and the defining equations among these
// new variables encode the fact that A * B = B * A = W.I
//
// We assume R = basering. The return value is a ring nR and in this ring
// is already declared an ideal moduli_eqns which is the quotient, 
// i.e. nR/moduli_eqns is the parametrising space.

proc moduliOfPotential(poly W, int r, int k)
{
	if( typeof(basering) != "ring" )
	{
		print("[moduliOfPotential] This is the wrong function for quotient rings. Exiting.");
		return();
	}
	
	if( W == 0 )
	{
		print("[moduliOfPotentials] Not a good idea to pass W = 0 in here yet.");
		return();
	}
	
	// Monomials of total weight <= k in the ring variables, as intvecs
	int numVars = nvars(basering);
	list monoms = monomialdict(k, numVars);
 	int Q = size(monoms);
 	
	// Get the coefficient vector of W and the monomials that occur
    poly xprod = 1;
    int i;
    for(i=1; i<=numVars; i++)
    {
        xprod = xprod * var(i);
    }
    matrix koffer = coef(W, xprod);
    
    // the first row of koffer contains the monomials and the second row the coeffs
    // store this as a list of monomial intvecs and the coeffs
    list mm_list, mc_list;
    for(i=1; i<=ncols(koffer); i++)
    {
    	mm_list = mm_list + list(leadexp(koffer[1,i])); // leadexp extracts intvec
    	mc_list = mc_list + list(number(koffer[2,i])); // number to force conversion
    }
	
    // The variables of our new ring are
    // f(z)(i)(j)(a) from f(1,2)(1..r)(1..r)(1..Q)
    // where z = 1 means top right block A, z = 2 means bottom left B,
    // i and j are the row and column resp., and a is the monomial index
    list newVarNames;
    int z,i,j,a;
    for( z=1; z<=2; z++ )
    {
    	for( i=1; i<=r; i++ )
    	{
    		for( j=1; j<=r; j++ )
    		{
    			for( a=1; a<=Q; a++ )
    			{
    				// WARNING: If you change the ordering here you should also change it in e.g.
    				// universalMF and other places where the order is used
    				string nv = "f(" + string(z) + ")(" + string(i) + ")(" + string(j) + ")(" + string(a) + ")";
    				newVarNames = newVarNames + list(nv);
    			}
    		}
    	}
    }
    //dbprint(printlevel, "[moduliOfPotential] Var names are f(z)(i)(j)(a) where 1 <= a <= " + string(Q) + " is an index into monoms");
    //dbprint(printlevel, "[moduliOfPotential] Number of moduli variables: " + string(size(newVarNames)));
    
    // Now create the new ring
    def RRR = basering;
    
	list nR_list;
	nR_list[1] = ringlist(RRR)[1]; // TOOD: minpoly
	nR_list[2] = newVarNames;
	
	intvec kk = (1..size(newVarNames));
    for(i=1; i<=size(newVarNames); i++)
    {
        kk[i] = 1;
    }
	nR_list[3] = list( list( "dp", kk ) );
	nR_list[4] = ideal(0);
	
	def nR = ring(nR_list);
    setring nR;
    
    // mm_list is already available
    list mc_list = fetch(RRR,mc_list);
    	
    // Find the ideal of moduli equations. The equations encode the
    // fact that A * B = W.I and B * A = W.I. Thus we have one equation
    // for each pair (i,j) of matrix positions and each monomial which is
    // either of total weight <= 2k or occurs in W
    
    // First make a list of such monomials (as intvecs)
    list eqn_monoms = monomialdict(2*k, numVars);
	
	// Go through the monomials in W and add them to eqn_monoms if they don't already occur
	for(i=1;i<=size(mm_list);i++)
	{
		if( sum(mm_list[i]) > 2*k )
		{
			eqn_monoms = eqn_monoms + list(mm_list[i]);
		}
	}
	
    list moduli_eqn_list;
    int u;
    
    for( i=1; i<=r; i++ )
    {
    	for( j=1; j<=r; j++ )
    	{
    		for( u=1; u<=size(eqn_monoms); u++ )
    		{
    			poly e1,e2;
    			intvec m = eqn_monoms[u];
    			
    			// At row i, col j, and for the monomial m the associated
    			// equation e1 asserts that the coeff of m in (A*B)_{i,j} equals the coeff
    			// of m in (W.I)_{i,j} and e2 asserts the same for B*A.
    		
    			// The coeff of m in (W.I)_{i,j} is given by cw, where
    			poly cw;
    			if( i != j )
    			{
    				cw = 0;
    			}
    			else
    			{
    				// On diagonal i = j
    				int a = indexinlist(mm_list,m);
    				
    				if( a == 0 )
    				{
    					// This monomial has zero coeff in W
    					cw = 0;
    				}
    				else
    				{
    					cw = mc_list[a];
    				}
    			}
    			
    			e1 = -cw;
    			e2 = -cw;
    			
    			// The coeff of m in (A*B)_{i,j} is a sum of numbers
    			// f(1)(i)(k)(c) * f(2)(k)(j)(d) where k ranges over 1...r
    			// and independently c,d range over indices between 1..Q
    			// with the condition that monoms[c] + monoms[d] = m.
    			// The coeff of m in (B*A)_{i,j} is a the same but with 1,2 interchanged
    			int k, ac, ad;
    			for( k=1; k<=r; k++ )
    			{
    				for( ac=1; ac<=Q; ac++ )
    				{
    					for( ad=1; ad<=Q; ad++ )
    					{
    						if( monoms[ac] + monoms[ad] == m )
    						{
    							e1 = e1 + f(1)(i)(k)(ac) * f(2)(k)(j)(ad);
    							e2 = e2 + f(2)(i)(k)(ac) * f(1)(k)(j)(ad);
    						}
    					}
    				}
    			}
    			     
    			moduli_eqn_list = moduli_eqn_list + list(e1,e2);
    		}
    	}
    }
    
    // convert to a row matrix
    matrix mat_eqn[1][size(moduli_eqn_list)];
    for( i=1; i<=size(moduli_eqn_list); i++ )
    {
    	mat_eqn[1,i] = moduli_eqn_list[i];
    }
    
    ideal moduli_eqns = ideal(mat_eqn);    
    export(moduli_eqns);
            
	setring RRR;
	return(nR);
}

////////////////////////////////////////////////////////////////////
// universalMF
//
// Given the same input data as moduliOfPotential returns the
// universal matrix factorisation of that moduli space. This is 
// defined over the ring R + nR where R is the basering and
// nR is the output of moduliOfPotential. 
//
// The return value is a ring uR and in this ring is already
// declared the ideal moduli_eqns from moduliOfPotential and
// a matrix factorisation U in supermatrix form (the universal MF).
// Note that this only satisfies U * U = W . I up to moduli_eqns.

proc universalMF(poly W, int r, int k)
{
	def RRR = basering;
	def nR = moduliOfPotential(W,r,k);
	def uR = RRR + nR; 
	
	setring(uR);
	
	// imap works to import stuff from RRR and nR to uR
	ideal moduli_eqns = imap(nR, moduli_eqns);
	poly W = imap(RRR,W);
	export(W);
	export(moduli_eqns);
	
	// Now construct the universal matrix U
	matrix A(1)[r][r]; matrix A(2)[r][r];
	
	list varNames = ringlist(RRR)[2];
 	int numVars = size(varNames);
 	
 	int i;
 	list monoms = monomialdict(k, numVars);
	int Q = size(monoms);

	int z,i,j,a;
    for( z=1; z<=2; z++ )
    {
    	for( i=1; i<=r; i++ )
    	{
    		for( j=1; j<=r; j++ )
    		{
    			for( a=1; a<=Q; a++ )
    			{
    				// The RRR variables come first, so if we call
    				// monomial with something of length <= numVars
    				// we get the monomial from RRR
    				A(z)[i,j] = A(z)[i,j] + f(z)(i)(j)(a) * monomial(monoms[a]);
    			}
    		}
    	}
    }
	
	matrix U = mfblocks(A(1),A(2));
	export(U);
		
	setring RRR;
	return(uR);
}

////////////////////////////////////////////////////////////////////
// moduliOfPotentialDet
//
// Given a polynomial ring R, a polynomial W in R, a rank r, and
// a degree cutoff k, returns a quotient ring which parametrises
// _determinantal_ matrix factorisations of W of rank r, with only
// monomials of total degree <= k.
//
// A determinantal MF is a pair (A,B) where det(A) = W and B is the 
// adjoint matrix of A.
//
// Details: Corresponding to each position (i,j) in A
// and monomial in R of total degree <= k we introduce a new variable f(i)(j)(a)
// where a is some index. These are the variables of our new ring, and the 
// defining equations among these new variables encode the fact that
// det(A) = W.
//
// We assume R = basering. The return value is a ring nR and in this ring
// is already declared an ideal moduli_eqns which is the quotient, 
// i.e. nR/moduli_eqns is the parametrising space.

proc moduliOfPotentialDet(poly W, int r, int k)
{
	if( typeof(basering) != "ring" )
	{
		print("[moduliOfPotentialDet] We cannot deal with quotient rings yet, sorry. Exiting.");
		return();
	}
	
	if( W == 0 )
	{
		print("[moduliOfPotentialDet] Not a good idea to pass W = 0 in here yet.");
		return();
	}
		 	
 	// Monomials of total weight <= k in the ring variables, as intvecs
 	int numVars = nvars(basering);
 	list monoms = monomialdict(k, numVars);
 	int Q = size(monoms);
	
    // Get the coefficient vector of W and the monomials that occur
    int i;
    poly xprod = 1;
    for(i=1; i<=numVars; i++)
    {
        xprod = xprod * var(i);
    }
    matrix koffer = coef(W, xprod);
    
    // the first row of koffer contains the monomials and the second row the coeffs
    // store this as a list of monomial intvecs and the coeffs
    list mm_list, mc_list;
    for(i=1; i<=ncols(koffer); i++)
    {
    	mm_list = mm_list + list(leadexp(koffer[1,i])); // leadexp extracts intvec
    	mc_list = mc_list + list(number(koffer[2,i])); // number to force conversion
    }
    
    // The variables of our new ring are
    // f(i)(j)(a) from f(1..r)(1..r)(1..Q)
    // i and j are the row and column resp., and a is the monomial index
    list newVarNames;
    int i,j,a;
    for( i=1; i<=r; i++ )
    {
    	for( j=1; j<=r; j++ )
    	{
    		for( a=1; a<=Q; a++ )
    		{
    			// WARNING: If you change the ordering here you should also change it in e.g.
    			// universalMFDet and other places where the order is used
    			string nv = "f(" + string(i) + ")(" + string(j) + ")(" + string(a) + ")";
    			newVarNames = newVarNames + list(nv);
    		}
    	}
    }
    
    // Now create the new ring
    def RRR = basering;
    
	list nR_list;
	nR_list[1] = ringlist(RRR)[1]; // TOOD: minpoly
	nR_list[2] = newVarNames;
	
	intvec kk = (1..size(newVarNames));
    for(i=1; i<=size(newVarNames); i++)
    {
        kk[i] = 1;
    }
	nR_list[3] = list( list( "dp", kk ) );
	nR_list[4] = ideal(0);
	
	def nR = ring(nR_list);
    setring nR;
    
    // mm_list is already available
    list mc_list = fetch(RRR,mc_list);
    
    // Find the ideal of moduli equations. The equations encode the
    // fact that det(A) = W. Thus we have one equation for each monomial
    // which occurs in det(A) or W. The former are all monomials of total
    // weight <= rk
    
    list monoms_rk = monomialdict(r*k,numVars);
 	int Q_rk = size(monoms_rk);
 	
 	// We need a list L of all intvecs v of length r whose entries
 	// are between 1 and Q, to index a choice of monomial of weight <= k
 	// for each row of A in the determinant calculation
 	list L = seqsfrominterval(Q, r);
 	
 	// Make a list L_sum which adds up all the exponents from L
 	// That is, each element L[i] chooses a sequence of r monomials of
 	// total weight <= k, namely monoms[L[i][1]], ..., monoms[L[i][r]]
 	// and L_sum[i] = monoms[L[i][1]] + ... + monoms[L[i][r]]
 	int j;
 	list L_sum;
 	intvec v;
 	for(i=1;i<=size(L);i++)
 	{
 		v = 0;
 		
 		for(j=1;j<=r;j++)
 		{
 			v = v + monoms[L[i][j]];
 		}
 		
 		L_sum[i] = v;
 	}
 	
    // To compute the determinant we first need all permutations of r
    list perms = SGroupintvecs(r);
    
    // Zero the list of coefficients of det(A)
    list det_coeff;
    for(i=1;i<=Q_rk;i++)
    {
    	det_coeff[i] = 0;
    }
    
    // Now for each permutation s = perms[i], and each sequence of choices 
    // of monomials L[j] (one per row) we get a contribution to the coefficient
    // of the monomial L_sum[j] in det(A)
    int i, j, u;
    intvec s;
    intvec mseq;
    intvec lsum;
    int b;
    poly t;
    
    for(i=1;i<=size(perms);i++)
    {
		s = perms[i];
	    
    	for(j=1;j<=size(L);j++)
    	{
    		mseq = L[j];
    		
    		t = (-1)^(LengthSymElement(s));
    		
    		// u is the row
	    	for(u=1;u<=r;u++)
    		{
    			t = t * f(u)(s[u])(mseq[u]);
	    	}	
	    	
	    	// Add t to the coefficient of L_sum[j] in det(A)
	    	b = indexinlist(monoms_rk, L_sum[j]);
	    	
	    	det_coeff[b] = det_coeff[b] + t;
	    }
    }
        
    // The monomials for which we need equations are monoms_rk plus the
    // monomials in W. Go through the monomials in W and add them
    // to eqn_monoms if they don't already occur
    list eqn_monoms = monoms_rk;
	
	for(i=1;i<=size(mm_list);i++)
	{
		if( sum(mm_list[i]) > r*k )
		{
			eqn_monoms = eqn_monoms + list(mm_list[i]);
		}
	}
	
	string deb;
	for(i=1;i<=size(monoms);i++)
	{
		deb = deb + "[" + string(monoms[i]) + "] ";
	}
	//dbprint(printlevel, "[moduliOfPotentialDet] monoms:" + deb);
	//dbprint(printlevel, "[moduliOfPotentialDet] det_coeff:" + string(det_coeff));

	// Now generate the moduli equations
	list moduli_eqn_list;
    int u;
    
	for( u=1; u<=size(eqn_monoms); u++ )
	{
		poly e;
    	intvec m = eqn_monoms[u];
    			
    	// At row i, col j, and for the monomial m the associated
    	// equation e asserts that the coeff of m in det(A) equals the coeff
    	// of m in W.
    		
    	// The coeff of m in W is given by cw, where
    	poly cw;
    	int a = indexinlist(mm_list,m);
    				
    	if( a == 0 )
    	{
    		// This monomial has zero coeff in W
    		cw = 0;
    	}
    	else
    	{
    		cw = mc_list[a];
    	}
    			
    	e = -cw;
    			
    	// The coeff of m in det(A) is similarly computed
    	int b = indexinlist(monoms_rk,m);
	    	
	    if( b != 0 )
	    {
	    	e = e + det_coeff[b];
		}
    			
    	moduli_eqn_list = moduli_eqn_list + list(e);
    }
    
    // convert to a row matrix
    matrix mat_eqn[1][size(moduli_eqn_list)];
    for( i=1; i<=size(moduli_eqn_list); i++ )
    {
    	mat_eqn[1,i] = moduli_eqn_list[i];
    }
    
    ideal moduli_eqns = ideal(mat_eqn);    
    export(moduli_eqns);
        
	setring RRR;
	return(nR);
}

////////////////////////////////////////////////////////////////////
// universalMFDet
//
// Given the same input data as moduliOfPotentialDet returns the
// universal matrix factorisation of that moduli space. This is 
// defined over the ring R + nR where R is the basering and
// nR is the output of moduliOfPotential. 
//
// The return value is a ring uR and in this ring is already
// declared the ideal moduli_eqns from moduliOfPotential and
// a matrix factorisation U in supermatrix form (the universal MF).
// Note that this only satisfies U * U = W . I up to moduli_eqns.

proc universalMFDet(poly W, int r, int k)
{
	def RRR = basering;
	def nR = moduliOfPotentialDet(W,r,k);
	def uR = RRR + nR; 
	
	setring(uR);
	
	// imap works to import stuff from RRR and nR to uR
	ideal moduli_eqns = imap(nR, moduli_eqns);
	poly W = imap(RRR,W);
	export(W);
	export(moduli_eqns);
	
	// Now construct the universal matrix U beginning with the upper
	// right block, which is what the moduli space parametrises
	matrix A[r][r]; 
	
 	int numVars = nvars(RRR);
 	list monoms = monomialdict(k,numVars);
	int Q = size(monoms);

	int i,j,a;
    for( i=1; i<=r; i++ )
    {
    	for( j=1; j<=r; j++ )
    	{
    		for( a=1; a<=Q; a++ )
    		{
    			// The RRR variables come first, so if we call
    			// monomial with something of length <= numVars
    			// we get the monomial from RRR
    			A[i,j] = A[i,j] + f(i)(j)(a) * monomial(monoms[a]);
    		}
    	}
    }
    
    // The other matrix B is the adjugate of A
	matrix B = adjoint(A);
	
	matrix U = mfblocks(A,B);
	export(U);
		
	setring RRR;
	return(uR);
}
	
////////////////////////////////////////////////////////////////////
// monomialdict
//
// Given integers k >= 1 and nv >= 1 returns all intvecs of length
// nv with sum <= k.

proc monomialdict(int k, int nv)
{
 	int i;
 	list monoms;
 	for(i=0; i<=k; i++)
 	{
		monoms = monoms + partitions(i, nv);
	}
	
	return(monoms);
}

////////////////////////////////////////////////////////////////////
// partitions
//
// Given integers N and M with N >= 0 and M > 0 returns all elements
// in NN^M which sum to N (here NN means non-negative integers) as a list
// of intvecs.

proc partitions(int N, int M)
{
	if( N < 0 || M <= 0 )
	{
		print("[partitions] Bad integers, exiting.");
		return();
	}
	
	if( M == 1 )
	{
		return(list(intvec(N)));
	}
	
	list P;
    int i, j;
    
    for( i=0; i<=N; i++)
    {
    	list miniP = partitions(N-i, M-1);
    
    	for( j = 1; j<=size(miniP); j++ )
    	{
    		P = P + list( concat_intvec( intvec(i), miniP[j] ) );
    	}
    }
    
    return(P);
}

////////////////////////////////////////////////////////////////////
// displayIndices

proc displayIndices(int k)
{
	int numVars = nvars(basering);
	list monoms = monomialdict(k, numVars);
 	
 	print("Monomials of total degree <= " + string(k));
	print("");
 	
 	int i;
 	for(i=1;i<=size(monoms);i++)
 	{
 		print(string(i) + "    " + string(monomial(monoms[i])) + "    [" + string(monoms[i]) + "]");
 	}
}

////////////////////////////////////////////////////////////////////
// seqsfrominterval
//
// Given an integer N >=1 and M >= 1 return the set of all sequences
// of integers of length M with entries in 1,...,N

proc seqsfrominterval(int N, int M)
{
	if( M == 1 )
	{
		list L;
		int i;
		for(i=1;i<=N;i++)
		{
			L = L + list(intvec(i));
		}
		return(L);
	}
	
	list L_small = seqsfrominterval(N, M-1);
	
	list L;
	int i,j,k;
	intvec v;
	
	for(i=1;i<=N;i++)
	{
		for(j=1;j<=size(L_small);j++)
		{
			v[1] = i;
			
			for(k=1;k<=M-1;k++)
			{
				v[k+1] = (L_small[j])[k];
			}
			
			L = L + list(v);
		}
	}
	
	return(L);
}

proc concat_intvec( intvec L, intvec R )
{
	int n = size(L);
	int m = size(R);
	intvec ret;
	
	int i;
	for(i=1;i<=n;i++)
	{
		ret[i] = L[i];
	}
	
	for(i=1;i<=m;i++ )
	{
		ret[n+i] = R[i];
	}
	
	return(ret);
}