blow.lib

version="1.0";
category="Miscellaneous";
info="
LIBRARY:   blow.lib  Pushforwards of matrix factorisations
AUTHOR:    Nils Carqueville, 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";
LIB "qmatrix.lib";
LIB "control.lib";

////////////////////////////////////////////////////////////////////
// USAGE GUIDE
//
// We expect our ambient ring to be of one of the types
//
//     0,(x(1..t),y(1..s)),dp    OR    (0,r),(x(1..t),y(1..s)),dp
// 
// for some integers t and s. The x-variables are the "internal" variables.
// We always expect "r" to be the param defining an extension ring.
//
// In the second case where we define a ring extension, in addition to defining
// the minpoly we expect that before using this library an additional variable
// "minblowpoly" is defined, which is the minimum polynomial of r over Q expressed
// in terms of the variable x(1). For example
//
//    ring rr=(0,r),(x(1..t),y(1..s)),dp;
//    minpoly = r^2 + 1;
//    poly minpolyblow = x(1)^2 + 1;
//
// If minpolyblow or r is undefined we proceed by ignoring the ring extension.
//
// We assume currently in dQ, deltaQ and related routines that the polynomial
// we are passed is a sum of all the x-variables to some specific power.
//
// NOTE: With Singular 4 the way minimum polynomials are handled has been changed,
// see http://www.singular.uni-kl.de/Manual/4-0-1/sing_2235.htm#SEC2313.

////////////////////////////////////////////////////////////////////
// GRADED MATRIX FACTORISATIONS
//
// IMPORTANT: Following KR we make use of the grading on our polynomial rings where each
// variable has degree 2. This is NOT imposed on the level of Singular code, it just affects
// our convention for grading vectors (which terminology will be explained presently). Given
// a polynomial p we write deg(p) for the usual degree, as understood by Singular, and |p| for
// the doubled grading, so |p| = 2 * deg(p). Given a graded module M, M{1} denotes the module
// with grading M{1}_i = M_{i-1}.
//
// 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.
//
// A matrix factorisation is represented by an odd supermatrix A giving the differential
// on some Z/2-graded free module X (over a polynomial ring). A _graded_ matrix factorisation
// is this differential together with an intvec g, whose length is equal to the number of
// columns in A, giving the "grading shift" on each free module in X. We require that the
// components a0 and a1 of A are then morphisms of graded modules of degree n+1 where the
// potential is of degree 2(n+1).
//
// For example, if X = X0 + X1 is a free module of rank 4, then g = (3 3 2 1) would express
// that as a graded module X0 = R{3} + R{3} and X1 = R{2} + R{1} where R is the polynomial
// ring with the doubled grading. Suppose that the potential W has |W| = 4, so n = 1. Then
// the differential A = (0, a1 newline a0, 0) would have to consist of a morphism
// a0: R{3}+R{3} -> R{2}+R{1} of degree n+1 = 2, so 
// |a0[1,1]| = 3, (i.e. this entry would have to be zero, because there are no odd degrees)
// |a0[1,2]| = 3,
// |a0[2,1]| = 4 (i.e. a polynomial f with deg(f) = 2)
// |a0[2,2]| = 4
//
// NOTE: All operations on graded MFs assume that A is square (i.e. both degrees have
// the same rank).
//
// NOTE: Bigradings
//
// In the case of bigradings (for HOMFLY homology) the first ring variable (call it "a")
// has bidegree (2,0) and the rest have bidegree (0,2). The grading shift on a free module of
// rank M is described by a list of size M (even if M = 1) whose entries are intvecs of size
// 2. So R{a,b} + R{c,d} is described by list(invec(a,b),intvec(c,d)).

//////////////////////////////////////////////////
// blowFLags
//
// Functions as a list of #define's.
// Change the return value, say of sanity_checks, to turn off sanity checks in all routines.
//
// NOTE: Singular's evaluation of a && b does not use shortcuts, it evaluates both a and
// b before doing the logical AND, so if one uses flags to avoid computation, one has 
// to nest ifs. 

proc blowFlags(string s)
{
    // Throughout the code, there are various checks on the integrity of data
    // (that certain things are morphisms, that MFs factorise the potential they
    // should, and so on). These can slow down the code dramatically, but help
    // locate bugs. Turn off for large computations.
    if( s == "sanity_checks" )
    {
        if( defined(blow_flag_sanity_checks) )
        {
            return(1); 
        }
        else
        {
            return(0); // default is 0
        }
    }
    
    // Turn on this flag if you want to optimise for small coefficients (of monomials)
    // rather than speed, in for example mfReduce.
    if( s == "small_coeffs" )
    {
        return(0); // default is 0
    }
    // DEBUG: As of Singular 4, the routine findConstantMin has a bug which
    // causes the interpreter to crash with a Segmentation fault. So until this
    // is fixed, DO NOT USE small_coeffs
    
    // Enables some checks only relevant if the calculations involved
    // are for link homology (e.g. maps in webCompileMorphism should
    // be degree +1)
    if( s == "link_homology_checks" )
    {
        return(0); // default is 0
    }
    
    if( s == "minbase_cutoff" )
    {
        return(5); // default is 5
    }
    
    // We make use of two storage options: ASCII and MPFILE. The flag use_mpfile
    // controls which one is used (some operating systems only support ASCII).
    // The default is to use ASCII since it is more compatible, but MPFILE is 
    // much faster.
    if( s == "use_mpfile" )
    {
        if( defined(blow_flag_use_mpfile) )
        {
            return(1);
        }
        else
        {
            return(0); // default is 0
        }
    }
    
    // Cutoff for how high to look for a power of ring variables belonging
    // to the Jacobi ideal
    if( s == "power_cutoff_jacobi" )
    {
        return(1000); // default is 1000
    }
    
    if( s == "factor_out_content" )
    {
        return(0); // default is 0
    }
    
    // We did not recognise this flag
    print("[blowFlags] Flag not recognised, returning zero.");
    
    return(0);
}

//////////////////////////////////////////////////////////////////////////////////
// numXVars - Returns the number of x variables (i.e. "nx")
//
// We assume that the ring is structured as described at the beginning of this document,
// so that the variables are x(1),..,x(t),y(1),...,y(s) for some integers t and s. This
// routine returns the integer t (possibly zero).
//////////////////////////////////////////////////////////////////////////////////

proc numXVars()
{
    list varNames = ringlist(basering)[2];
    
    int num = 0;
    int i;
    for(i=1;i<=size(varNames);i++)
    {
        if( varNames[i][1] == "x" )
        {
            num++;
        }
    }
    
    return(num);
}

//////////////////////////////////////////////////////////////////////////////////
// numYVars - Returns the number of y variables (i.e. "ny")
//
// We assume that the ring is structured as described at the beginning of this document,
// so that the variables are x(1),..,x(t),y(1),...,y(s) for some integers t and s. This
// routine returns the integer s (possibly zero).
//////////////////////////////////////////////////////////////////////////////////

proc numYVars()
{
    list varNames = ringlist(basering)[2];
    
    int num = 0;
    int i;
    for(i=1;i<=size(varNames);i++)
    {
        if( varNames[i][1] == "y" )
        {
            num++;
        }
    }
    
    return(num);
}

//////////////////////////////////////////////////////////////////////////////////
// "ringwithoutyvars" returns the ambient ring with y-variables removed
// NOTE: We just return the ring, we do not change the current basering
// NOTE: This assumes our only ring variables are xs and ys
// NOTE: This does get called by linkhom.lib with numY = 0, i.e. nothing to do.
//////////////////////////////////////////////////////////////////////////////////

proc ringWithoutYVars()
{
    def RRR = basering;
    int numX = numXVars();
    int numY = numYVars();
    
    if( numX == 0 )
    {
        print("[ringWithoutYVars] Cannot create a ring with no variables, failing.");
        return();
    }
    
    // Define new ring -- exactly the same as RRR, but without the y-variables:
    // ringlist(RRR)[1] will be the characteristic followed by (possibly) the names of params.
    if(size(ringlist(RRR)[1]) > 1)
    {
        int i,j,k;
        list L = ringlist(RRR);
        list newVars;
        for(i=1;i<=numX;i++)
        {
            string nv = "x(" + string(i) + ")";
            newVars = newVars + list(nv);
            kill nv;
        }
        L[2] = newVars;
        kill newVars;
        
        intvec kk = (1..numX);

        for(i=1; i<=numX; i++)
        {
            kk[i] = 1;
        }
        L[3][1][2] = kk;
        L[1][4][1] = 0;     // set minpoly to zero, put it back to the right value below...
        def nR = ring(L);

        // Now complete the correct definition of nR by specifying the right minpoly: 
        if( defined(minpolyblow) && defined(r) )
        {
            setring nR;
            poly minpolyblow = fetch(RRR,minpolyblow); // Use fetch so the minpoly stays using var(1)
            poly z = subst(minpolyblow,var(1),r);
            number nu = leadcoef(z);
            minpoly = nu;
            poly minpolyblow = fetch(RRR,minpolyblow);
	        export(minpolyblow);
    	    // NOTE: The double definition of minpolyblow is because Singular nukes
        	// all variables in our ring when we set minpoly
        }
        else
        {
            print("[ringWithoutYVars] Ring param other than r, or unspecified minpolyblow. Please see the Usage Guide.");
        }
    }
    else
    {
        ring nR = char(RRR), x(1..numX), dp;
    }

    def newRing = basering;
    
    // Before returning the new ring, set the ambient ring as current:
    setring RRR;
    
    return(newRing);
}


//////////////////////////////////////////////////////////////////////////////////
// "ringWithoutXVars" returns the ambient ring with x-variables removed
// NOTE: We just return the ring, we do not change the current basering
//////////////////////////////////////////////////////////////////////////////////

proc ringWithoutXVars()
{
    def RRR = basering;
    int numX = numXVars();
    int numY = numYVars();
    
    if( numY == 0 )
    {
        print("[ringWithoutXVars] Cannot create a ring with no variables, failing.");
        return();
    }
    
    // Define new ring -- exactly the same as RRR, but without the x-variables:
    // ringlist(RRR)[1] will be the characteristic followed by (possibly) the names of params.
    if(size(ringlist(RRR)[1]) > 1)
    {
        int i,j,k;
        list L = ringlist(RRR);
        
        list newVars;
        for(i=1;i<=numY;i++)
        {
            string nv = "y(" + string(i) + ")";
            newVars = newVars + list(nv);
            kill nv;
        }
        L[2] = newVars;
        kill newVars;
        
        intvec kk = (1..numY);

        for(i=1; i<=numY; i++)
        {
            kk[i] = 1;
        }
        L[3][1][2] = kk;
        L[1][4][1] = 0;     // set minpoly to zero, put it back to the right value below...
        def nR = ring(L);

        // Now complete the correct definition of nR by specifying the right minpoly: 
        if( defined(minpolyblow) && defined(r) )
        {
            setring nR;
            poly minpolyblow = fetch(RRR,minpolyblow);
            poly z = subst(minpolyblow,var(1),r);
            number nu = leadcoef(z);
            minpoly = nu;
            poly minpolyblow = fetch(RRR,minpolyblow);
	        export(minpolyblow);
    	    // NOTE: The double definition of minpolyblow is because Singular nukes
        	// all variables in our ring when we set minpoly
        }
        else
        {
            print("[ringWithoutXVars] Ring param other than r, or unspecified minpolyblow. Please see the Usage Guide.");
        }
    }
    else
    {
        ring nR = char(RRR), y(1..numY), dp;
    }    
    
    def newRing = basering;
    
    // Before returning the new ring, set the ambient ring as current:
    setring RRR;
    
    return(newRing);
}


//////////////////////////////////////////////////////////////////////////////////
// "dimAlgebraOverInternalVariables" Returns the dimension of the quotient of the
// polynomial ring in the internal variables (as returned by ringWithoutYVars) by the
// specified ideal (given as an ideal in the basering).
//////////////////////////////////////////////////////////////////////////////////

proc dimAlgebraOverInternalVariables(ideal J)
{	
    def RRR = basering;
    def nR = ringWithoutYVars();
    setring nR;
	
    ideal J = imap(RRR,J);
    ideal Jstd = std(J);
    ideal ba = kbase(Jstd);
    int n = size(ba);

    setring RRR;

    return(n);
}

//////////////////////////////////////////////////////////////////////////////////
// "deltaInflationList" Returns a list containing, in order, the matrix representing
// the operator [g_i,x(i)^b] for b = 1,...,2*exponent - 1.
//////////////////////////////////////////////////////////////////////////////////

proc deltaInflationList(int i, int exponent, ideal J)
{
    int prl = printlevel;
    
    // Compute a k-basis of J:
    def RRR = basering;
    def nR = ringWithoutYVars();
    setring nR;
    
    ideal J = imap(RRR,J);
    ideal Jstd = std(J);
    ideal basis = kbase(Jstd);
    module BB = reduce(basis,Jstd);

    int jdim = size(basis);
    int b, k;
    list deltaActions;
    
    for(b=1; b<= 2*exponent - 1; b++)
    {
        // Compute the matrix representing the operator [g_i, x(i)^b] on I = nR/J:
        list L;
        for(k=1; k<=jdim; k++)
        {
            // First compute g_i( x(i)^b * - ) on our basis element:
            poly f1 = ( x(i)^b * basis[k] - reduce(x(i)^b * basis[k],std(x(i)^exponent)) )/(x(i)^exponent);

            // Next compute x(i)^b * g_i(-):
            poly f2 = x(i)^b * ( basis[k] - reduce(basis[k],std(x(i)^exponent)) )/(x(i)^exponent);

            // To our output matrix we add the expression of the commutator f1 - f2 in the
            // chosen basis of the algebra I:
            module ff = reduce(f1 - f2,Jstd);
            matrix MM = matrix(reduce(lift(BB,ff),std(syz(BB))));
            L = L + list(MM);
            kill f1,f2,ff,MM;
        }

        matrix U = L[1];
        for(k=2; k<=jdim; k++)
        {
            U=concat(U,L[k]);
        }
                
        deltaActions = deltaActions + list(U);
        kill U, L;
    }
    
    setring RRR;
    list deltaActions = imap(nR,deltaActions);
    
    return(deltaActions);
}

//////////////////////////////////////////////////////////////////////////////////
// "variableInflationList" Returns a list containing, in order, the matrix representing
// the action of each of the internal variables on k[x]/J (the internal algebra).
//////////////////////////////////////////////////////////////////////////////////

proc variableInflationList(ideal J)
{
    int prl = printlevel;
    
    def RRR = basering;
    def nR = ringWithoutYVars();
    setring nR;
    
    ideal Jx = std(imap(RRR,J));
    ideal basis = kbase(Jx);
    int jdim = size(basis);
    
    if( jdim == 0 )
    {
        print("[variableInflationList] Internal algebra is zero-dimensional, exiting.");
        return();
    }
    
    module BB = reduce(basis,Jx);
    module syzBB = std(syz(BB)); 
    
    int i, k;
    list variableActions;
    int numX = nvars(nR);
    
    for(i=1; i<=numX; i++)
    {
        list L;
        for(k=1; k<=jdim; k++)
        {
            module ff = reduce(x(i)*basis[k],Jx);
            matrix MM = matrix(reduce(lift(BB,ff), syzBB));
            L = L + list(MM);
        }
        
        matrix U = L[1];
        for(k=2; k<=jdim; k++)
        {
            U=concat(U,L[k]);
        }
        
        variableActions = variableActions + list(U);
        kill U, L;
    }
    
    setring RRR;
    list variableActions = imap(nR,variableActions);
    
    return(variableActions);
}

//////////////////////////////////////////////////////////////////////////////////
// "poly2matrix" blows up a polynomial into a matrix: x-monomials are replaced by
// the matrices that represent the former's action on the algebra k[x-variables]/J.
//
// NOTE: We assume that the algebra k[x-variables]/J is finite-dimensional.
//
// We must be passed a list "variableActions" which gives, in order, the matrices
// representing the action of multiplication by each internal variable on the
// internal algebra (as given by variableInflationList).
//
//////////////////////////////////////////////////////////////////////////////////

proc poly2matrix(poly tepo, list variableActions)
{
    // First compute the x-monomials and their coefficients in tepo:
    int i;
    poly xprod = 1;
    int numX = numXVars();
    
    for(i=1; i<=numX; i++)
    {
        xprod = xprod * x(i);
    }
    matrix koffer = coef(tepo, xprod);

    int u,k;
    int jdim = ncols(variableActions[1]);
    matrix inflation[jdim][jdim];
    
    for(u=1;u<=ncols(koffer);u++)
    {
        // Each column of koffer represents an x-monomial, let us now factorise this.
        // The list powers contains, in order, the power of each internal variable
        // in monom + 1 (since we multiply by xprod to make sure all vars appear).
        poly monom = koffer[1,u];
        intvec powers = factorize(monom*xprod,2)[2];
        matrix monomAction = unitmat(jdim);
        
        for(k=1;k<=numX;k++)
        {
            monomAction = monomAction * power(variableActions[k],powers[k]-1);
        }
        
        inflation = inflation + koffer[2,u] * monomAction;
        
        kill monom, powers, monomAction;
    }
    
    return(inflation);
}

////////////////////////////////////////////////////////////////////
// mfInflate
//
// We are given a matrix A interpreted as the differential on a matrix factorisation X
// of a potential W in our ring, together with a ring variable intvar and a power N.
// We assume that W is a polynomial in the set of variables with intvar removed, and
// we return the differential on the tensor product X otimes Q[intvar]/(intvar^N).

proc mfInflate(matrix A, poly intvar, int N)
{
    def RRR = basering;
    
    // Convert everything into a more suitable ring. We go through the original
    // list of ring variables, rename everything apart from intvar to y(i)'s (in
    // ascending order) and rename intvar to x(1).
    
    list rlist = ringlist(RRR);
    list varlist = rlist[2];
    
    // Create the new list of variables, where intvar is renamed to x(1)
    // and the other variables are enumerated as y(i)'s.
    list newvar;
    int i;
    int ycount;
    for(i=1;i<=size(varlist);i++)
    {
        if( string(intvar) == varlist[i] )
        {
            newvar = newvar + list("x(1)");
        }
        else
        {
            ycount++;
            string s = "y(" + string(ycount) + ")";
            newvar = newvar + list(s);
            kill s;
        }
    }
    
    list newringList;
    newringList[1] = rlist[1];
    newringList[2] = newvar;
    newringList[3] = rlist[3];
    newringList[4] = rlist[4];
    kill rlist;
    
    // Fix the variable weighting:
    intvec kk = (1..size(newvar));
    for(i=1; i<=size(newvar); i++)
    {
        kk[i] = 1;
    }
    newringList[3][1][2] = kk;
    kill newvar, kk;
    
    // Create our new ring:
    def nR = ring(newringList);
    kill newringList;
    setring nR;
    
    // Now we can call mablow:
    matrix A = fetch(RRR,A);
    ideal J = x(1)^N;
    matrix Ablow = mablow(A,J);
    
    setring RRR;
    matrix final = fetch(nR,Ablow);
        
    return(final);
}

//////////////////////////////////////////////////////////////////////////////////
// "mablow" blows up a matrix by blowing up all its entries using poly2matrix.
// The matrix M may be non-square.
//////////////////////////////////////////////////////////////////////////////////

proc mablow(matrix M, ideal J)
{
    int benchmark = 1; // Set to 1 to output benchmarking data

    // Compute dimension of the algebra k[x-vars]/J:
    int n = dimAlgebraOverInternalVariables(J);

    // Roughly speaking mablow takes "difficulty" ms to complete. We only give
    // debugging output if this is likely to take > 10 seconds.
    int difficulty = (nrows(M) * ncols(M) * n) div 3;
    if( difficulty > 10000 )
    {
        if( benchmark )
        {
            system("--ticks-per-sec",1000);
            int timeElapsed = timer;
        }
        
        dbprint(printlevel, "[mablow] Inflating " + string(nrows(M)) + "x" + string(ncols(M)) + " matrix by " + string(n) + "-dim algebra");
    }

    // Compute the inflation of individual variables:
    list variableActions = variableInflationList(J);
    
    // Define L to be an appropriately indexed list of blown-up matrices:
    int i1,j1,i2,j2,i;
    
    int ncolsM = ncols(M);
    int nrowsM = nrows(M);
    
    list e,L;
    for(i=1; i<=nrowsM; i++)
    {
        L[i] = e;
    }
    for(i1=1; i1<=nrowsM; i1++) // row
    {
        for(j1=1; j1<=ncolsM; j1++) // column
        {
            matrix PM = poly2matrix(M[i1,j1], variableActions);
            L[i1][j1] = PM;
            kill PM;
        }
    }

    matrix A[nrowsM*n][ncolsM*n];

    for(i1=1; i1<=nrowsM; i1++)
    {
        for(j1=1; j1<=n; j1++)
        {
            for(i2=1; i2<=ncolsM; i2++)
            {
                for(j2=1; j2<=n; j2++)
                {
                    A[(i1-1)*n + j1, (i2-1)*n + j2] = L[i1][i2][j1,j2];
                }
            }
        }
    }
    
    if( benchmark && difficulty > 10000 )
    {
        timeElapsed = timer - timeElapsed;
        dbprint(printlevel, "[mablow]            elapsed time " + string(timeElapsed) + "ms.");
    }

    return(A);
}

//////////////////////////////////////////////////////////////////////////////////
// "mablowGrading" returns the grading vector for an inflated MF, given the 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. Obviously this only makes
// sense if J is a graded ideal.
//////////////////////////////////////////////////////////////////////////////////

proc mablowGrading(intvec g, ideal J, int N)
{
    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);
        
    intvec r;   
    
    // There is an overall grading shift given by the sum of the degrees of
    // the generators of J, plus an additional N + 1 for each generator. 
    int k;
    int E;
    for(k=1;k <= size(J);k++)
    {
        E = E - 2 * deg(J[k]) + N + 1;
    }
    
    int i,j;
    for(i=1;i <= size(g); i++)
    {
        for(j=1;j<=jdim;j++)
        {
            r[(i-1)*jdim + j] = g[i] + 2 * deg(BB[j]) + E; 
        }
    }
    
    setring RRR;
    
    return(r);
}

//////////////////////////////////////////////////////////////////////////////////
// "poly2matrix_delta"
//
// Let us write nR for the basering without its y-variables.
//
// Given a polynomial in the x- and y-variables and an integer i, we return a matrix
// in only the y-variables, where the x-monomials are "inflated" by their action on
// the algebra I = nR/J.
//////////////////////////////////////////////////////////////////////////////////

proc poly2matrix_delta(poly tepo, int i, int exponent, list variableActions, list deltaActions)
{
    // First compute the x-monomials and their coefficients in tepo:
    int t;
    poly xprod = 1; 
    int numX = numXVars();
    for(t=1; t<=numX; t++)
    {
        xprod = xprod * x(t);
    }
    matrix koffer = coef(tepo, xprod);

    int u,k;
    int jdim = ncols(variableActions[1]);
    matrix inflation[jdim][jdim];
    
    for(u=1;u<=ncols(koffer);u++)
    {
        // Each column of koffer represents an x-monomial, let us now factorise this.
        // The list powers contains, in order, the power of each internal variable
        // in monom + 1 (since we multiply by xprod to make sure all vars appear).
        poly monom = koffer[1,u];
        intvec powers = factorize(monom*xprod,2)[2];
        matrix monomAction = unitmat(jdim);
        
        for(k=1;k<=numX;k++)
        {
            // If k is not equal to i, then we just take the usual inflation,
            // but for k = i we do the "delta" inflation:
            if( k != i )
            {
                monomAction = monomAction * power(variableActions[k],powers[k]-1);
            }
            else
            {
                // We have to multiply by a factor corresponding to x(i)^(powers[k]-1),
                // which is the action on the internal algebra of the operator
                // [g_i,x(i)^(powers[k]-1)], as recorded in deltaActions. If the exponent
                // is between 1 and 2 * exponent - 1 we use the array. Otherwise the
                // operator is zero, so multiply by zero.

                if( powers[k] - 1 >= 1 && powers[k] - 1 <= 2*exponent - 1)
                {
                    monomAction = monomAction * deltaActions[powers[k]-1];
                }
                else
                {
                    monomAction = 0 * monomAction;
                }
            }
        }
        
        inflation = inflation + koffer[2,u] * monomAction;
        
        kill monom, powers, monomAction;
    }
    
    return(inflation);
}

/////////////////////////////////////////////////////////////////////////////////////////
// "mablow_delta" blows up a matrix by blowing up all its entries using poly2matrix_delta.
//
// exponent is the integer such that g_i is division without remainder by x(i)^exponent.
/////////////////////////////////////////////////////////////////////////////////////////

proc mablow_delta(matrix M, int i, int exponent, ideal J)
{
    // Compute dimension of the algebra k[x-vars]/J:
    int n = dimAlgebraOverInternalVariables(J);
 
    // Compute the ordinary inflation of individual variables:
    list variableActions = variableInflationList(J);
    
    // Compute the delta inflation [g_i,x(i)^b] for all integers b:
    list deltaActions = deltaInflationList(i,exponent,J);
 
    // Define L to be an appropriately indexed list of blown-up matrices:
    int i1,j1,i2,j2,t;
    int s = ncols(M);
 
    list e,L;
    for(t=1; t<=s; t++)
    {
        L[t] = e;
    }
    for(i1=1; i1<=s; i1++)
    {
        for(j1=1; j1<=s; j1++)
        {
            L[i1] = insert(L[i1], poly2matrix_delta(M[i1,j1], i, exponent, variableActions, deltaActions),j1-1);
        }
    }   
    
    matrix A[s*n][s*n];

    for(i1=1; i1<=s; i1++)
    {
        for(j1=1; j1<=n; j1++)
        {
            for(i2=1; i2<=s; i2++)
            {
                for(j2=1; j2<=n; j2++)
                {
                    A[(i1-1)*n + j1, (i2-1)*n + j2] = L[i1][i2][j1,j2];
                }
            }
        }
    }

    return(A);
}

/////////////////////////////////////////////////////////////////////////////////////
// "SGroupintvecs" produces the n! permutations of (1,...,n) as a list of intvecs.
/////////////////////////////////////////////////////////////////////////////////////

proc SGroupintvecs(int n)
{
    list L;
    int i;
    
    intmat M = SymGroup(n);
    
    for(i=1; i<=nrows(M); i++)
    {
        intvec c = M[i,1..n];
        L[i] = c;
    }

    return(L);
}


/////////////////////////////////////////////////////////////////////////////////////
// "deltaQ" is delta Q^{wedge n} (where below we have F=Q).
//
// q is an intvec containing the exponents responsible for defining the operators g_i,
// so that g_i is division by x(i)^q[i] without remainder.
/////////////////////////////////////////////////////////////////////////////////////

proc deltaQ(matrix F, intvec q)
{
    //LIB "matrix.lib";
    
    // The ideal whose quotient gives the algebra I is generated by the internal
    // variables to the powers given by q:
    int i;
    ideal J;
    int numX = numXVars();
    
    for( i = 1; i <= numX; i++ )
    {
        J[i] = x(i)^(q[i]);
    }

    int n = dimAlgebraOverInternalVariables(J);
    int prl = printlevel;

    // The blown up matrix will be of dimension n x ncols(f):
    int gg = ncols(F) * n;

    list S = SGroupintvecs(numX);

    int j;    
    matrix Mi[gg][gg];

    // Precompute all the deltas:
    list deltas; 
    for(j = 1; j <= numX; j++)
    {
        deltas[j] = mablow_delta(F, j, q[j], J);
    }
    
    for(i=1; i<=size(S); i++)
    {
        intvec perm = S[i];
        matrix Ma = (-1)^(LengthSymElement(perm)) * unitmat(gg);
        
        for(j=1; j<=numX; j++)
        {
            Ma = Ma * deltas[perm[j]];
        }
        
        Mi = Mi + Ma; 
    }

    // Multiply by the scaling factor:
    number ffac = 1;
    for(i=1; i<=numX; i++)
    {
        ffac = ffac / i;
    }
    Mi = ffac * Mi;

    return(Mi);
}


/////////////////////////////////////////////////////////////////////////////////////
// "dQ" is the symmetrised product of null homotopies in the expression
// for the idempotent. The matrix f is the differential of the matrix factorisation
// and q, C refer to the reduction to null-homotopies with respect to powers of
// variables, see "Computing Khovanov-Rozansky homology and defect fusion".
/////////////////////////////////////////////////////////////////////////////////////

proc dQ(matrix f, intvec q, matrix C)
{
    //LIB "matrix.lib";
    poly detC = det(C);
    
    // The ideal whose quotient gives the algebra I is generated by the internal
    // variables to the powers given by q:
    int i;
    ideal J;
    int numX = numXVars();
    
    for( i = 1; i <= numX; i++ )
    {
        J[i] = x(i)^(q[i]);
    }

    int n = dimAlgebraOverInternalVariables(J);
    
    // The blown up matrix will be of dimension n x ncols(f):
    int gg = ncols(f) * n;

    list S = SGroupintvecs(numX);
    int j;

    matrix Mi[gg][gg];

    // Precompute all the mablows:
    list dels; 
    for(j = 1; j <= numX; j++)
    {
        dels[j] = mablow(diff(f,x(j)), J);
    }

    for(i=1; i<=size(S); i++)
    {
        intvec perm = S[i];
        matrix Ma = (-1)^(LengthSymElement(perm)) * unitmat(gg);    
        
        for(j=1; j<=numX; j++)
        {
            Ma = Ma * dels[perm[j]];
        }
        
        Mi = Mi + Ma; 
    }

    // Multiply with detC:
    Mi = mablow( detC * unitmat(ncols(f)), J) * Mi;

    // Multiply by the scaling factor:
    number ffac = 1;
    for(i=1; i<=numX; i++)
    {
        ffac = ffac / i;
    }
    Mi = ffac * Mi;
    
    return(Mi);
}


/////////////////////////////////////////////////////////////////////////////////////
// "varTrans" gives a list L with two entries:
// -- L[1]: an intvec q used to define g_i as division by x(i)^q[i] without remainder
// -- L[2]: a matrix C such that x(i)^q[i] = sum_j C_{i,j} partial_j W
// This is possible if and only if the Jacobi algebra of W is finite-
// dimensional, and we assume this is the case (otherwise this proc
// may not terminate).
/////////////////////////////////////////////////////////////////////////////////////

proc varTrans(poly W)
{
    def RRR = basering;
    def nR = ringWithoutYVars();
    setring nR;

    poly W = imap(RRR, W);
    ideal jaci = jacob(W);
    ideal J = std(jaci);
    int numVars = nvars(nR);

    int i,j,expo;
    intvec qvec;
    poly K;

    // Find the smallest possible integers qvec[i] such that x(i)^qvec[i] is in jacob(W):
    for(i=1; i<=numVars; i++)
    {
        int expo = 1;
        while( reduce( x(i)^expo, J ) != 0 )
        {
            expo++;
        }
        
        qvec[i] = expo;
    }

    // Define transformation matrix C:
    matrix C[numVars][numVars];
    for(i=1; i<=numVars; i++)
    {
        for(j=1; j<=numVars; j++)
        {
            C[j,i] = lift( jaci, x(j)^(qvec[j]) )[i,1];
        }
    }
    
    // Check that C does transform the partials to powers of the variables
    for(i=1;i<=numVars;i++)
    {
        poly f;
        int j;
        for(j=1;j<=numVars;j++)
        {
            f = f + C[i,j] * diff(W,x(j));
        }

        if( f != x(i)^(qvec[i]) )
        {
            print("[varTrans] Failed to compute transformation matrix correctly, exiting.");
            return();
        }
    }

    setring RRR;
    matrix C = imap(nR, C);
    list L = qvec, C;
    
    return(L);
}


/////////////////////////////////////////////////////////////////////////////////////
// "constantpart(P)" computes the constant part of a polynomial P.
/////////////////////////////////////////////////////////////////////////////////////

proc constantpart(poly P)
{
    int i;
    int n = nvars(basering);
    matrix C;
    C[1,1] = P;
    for(i=1; i<=n; i++)
    {
        C = coeffs( C[1,1], var(i) );
    }
    return(C[1,1]);
}


/////////////////////////////////////////////////////////////////////////////////////
// "permuterows(M, a, b)" permutes the rows a and b of a matrix M.
/////////////////////////////////////////////////////////////////////////////////////

proc permuterows(matrix M, int a,b)
{
    int i;
    int s = ncols(M);
    for(i=1; i<=s; i++)
    {
        poly pa = M[a,i];
        poly pb = M[b,i];
        M[a,i] = pb;
        M[b,i] = pa;
    }
    return(M);
}


/////////////////////////////////////////////////////////////////////////////////////
// "permutecols(M, a, b)" permutes the columns a and b of a matrix M.
/////////////////////////////////////////////////////////////////////////////////////

proc permutecols(matrix M, int a,b)
{
    int i;
    int s = nrows(M);
    for(i=1; i<=s; i++)
    {
        poly pa = M[i,a];
        poly pb = M[i,b];
        M[i,a] = pb;
        M[i,b] = pa;
    }
    return(M);
}

// blockmat(A,B,C,D) gives the block-matrix {{A,B},{C,D}} if the square matrices A,B,C,D all have the same size.

proc blockmat(matrix block1,matrix block2,matrix block3,matrix block4)
{
    int i,j;
    int r=ncols(block1);
    matrix m[2*r][2*r];
    for(i=1;i<=r;i++){for(j=1;j<=r;j++){m[i,j]=block1[i,j];}}
    for(i=1;i<=r;i++){for(j=r+1;j<=2*r;j++){m[i,j]=block2[i,j-r];}}
    for(i=r+1;i<=2*r;i++){for(j=1;j<=r;j++){m[i,j]=block3[i-r,j];}}
    for(i=r+1;i<=2*r;i++){for(j=r+1;j<=2*r;j++){m[i,j]=block4[i-r,j-r];}}
    return(m);
}

////////////////////////////////////////////////////////////////////
// mfblocks
//
// Given square matrices A, B of the same size creates a block matrix {{0,A},{B,0}}

proc mfblocks(matrix block1, matrix block2)
{
	if( nrows(block1) != ncols(block1) || nrows(block2) != ncols(block2) )
	{
		print("[mfblocks] Blocks are not square, exiting.");
        return();
    }
    
    if( nrows(block1) != nrows(block2) )
    {
    	printf("[mfblocks] Blocks are not the same size, exiting.");
    }
        
    int n = nrows(block1);
    
    matrix z[n][n];
	matrix Y = blockmat(z,block1,block2,z);
	return(Y);
}


// Modification of blockmat to accommodate nonsquare entries (consider replacing blockmat with this).
// We assume only that a block matrix makes sense, i.e. that in the order given the matrices fit
// together into a matrix. Here is the ASCII diagram:
//
//       a        b
//   ----------|-------|
//   |         |       |
//   |    1    |   2   |
// c |         |       |
//   |---------|-------|
// d |     3   |   4   |
//   |---------|-------|
//
// where we allow a, b, c, d to be arbitrary

proc blockmat_nsq(matrix block1,matrix block2,matrix block3,matrix block4)
{
    int i,j;
    int a = ncols(block1);
    int b = ncols(block2);
    int c = nrows(block1);
    int d = nrows(block3);
    
    // Sanity check
    if( a != ncols(block3) || b != ncols(block4) || c != nrows(block2) || d != nrows(block4) )
    {
        print("[blockmat] Incorrect block sizes, exiting.");
        return();
    }
    
    int nc=ncols(block1);
    int nr=nrows(block1);
    matrix m[c + d][a + b];

    for(i=1;i<=c;i++){for(j=1;j<=a;j++){m[i,j]=block1[i,j];}}
    for(i=1;i<=c;i++){for(j=a+1;j<=a+b;j++){m[i,j]=block2[i,j-a];}}
    for(i=c+1;i<=c+d;i++){for(j=1;j<=a;j++){m[i,j]=block3[i-c,j];}}
    for(i=c+1;i<=c+d;i++){for(j=a+1;j<=a+b;j++){m[i,j]=block4[i-c,j-a];}}
    return(m);
}


// To get blocks from block square matrix returns list (1-4 of blocks):

proc blocksq(matrix m)
{
    int n=nrows(m);
    int intnon2=n div 2;
    number non2=number(n) / 2;
    if(n!=ncols(m))
    {
        printf("Matrix should be square");
        return();
    }
    if(non2!=intnon2)
    {
        printf("Matrix should have an even number of rows/columns");
        return();
    }
    matrix block1[intnon2][intnon2]=submat(m,1..intnon2,1..intnon2);
    matrix block2[intnon2][intnon2]=submat(m,1..intnon2,intnon2+1..n);;
    matrix block3[intnon2][intnon2]=submat(m,intnon2+1..n,1..intnon2);;
    matrix block4[intnon2][intnon2]=submat(m,intnon2+1..n,intnon2+1..n);;
    list L=block1,block2,block3,block4;
    kill n,intnon2,non2,block1,block2,block3,block4;
    return(L);
}

////////////////////////////////////////////////////////////////////
// zeromat - Returns a square zero matrix of the specified size if
// called with one param, and a zero matrix of size n x m if called
// with two params n, m.
////////////////////////////////////////////////////////////////////

proc zeromat(int n, list #)
{
    if( size(#) == 0 )
    {
        matrix M[n][n];
        return(M);
    }
    else
    {
        int m = #[1];
        matrix M[n][m];
        return(M);
    }
}

////////////////////////////////////////////////////////////////////
// mfSuspend
// 
// Returns the suspension of the given matrix factorisation. If called
// with no additional params we suspend once, otherwise we suspend
// the number of times indicated by the second param.

proc mfSuspend(matrix A, list #)
{
    int power = 1;
    
    if( size(#) != 0 ){ power = #[1]; }
    
    list blocks = extractblockmat(A);
    matrix a0 = blocks[3]; matrix a1 = blocks[2];
    
    matrix z[ncols(A) div 2][ncols(A) div 2];
    
    if( power mod 2 == 0 )
    {
        // The power is even, so return the original MF:
        return(A);
    }
    else
    {
        // The power is odd, so suspend:
        matrix Asusp = blockmat(z, -1 * a0, -1 * a1, z);
        return(Asusp);
    }
}

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

////////////////////////////////////////////////////////////////////
// mfSuspendMorph
// 
// Returns the suspension of the given even morphism of matrix factorisations. If called
// with no additional params we suspend once, otherwise we suspend
// the number of times indicated by the second param.

proc mfSuspendMorph(matrix A, list #)
{
    int power = 1;
    
    if( size(#) != 0 ){ power = #[1]; }
    
    list blocks = extractblockmat(A);
    matrix a0 = blocks[1]; matrix a1 = blocks[4];
    
    matrix z[nrows(A) div 2][ncols(A) div 2];
    
    if( power mod 2 == 0 )
    {
        // The power is even, so return the original MF:
        return(A);
    }
    else
    {
        // The power is odd, so suspend:
        matrix Asusp = blockmat_nsq(a1, z, z, a0);
        return(Asusp);
    }
}

////////////////////////////////////////////////////////////////////
// mfStrictifyIdempotent
//
// Takes as input an odd block matrix A interpreted as the differential on a Z/2-graded free module X,
// and an even block matrix U interpreted as an endomorphism of X, satisfying U * U = U up to homotopy.
// We return, where possible, an even block matrix E representing an endomorphism of X homotopy
// equivalent to U which is strictly idempotent, meaning that E * E = E as chain maps.
//
// Here "where possible" means when (U * U - U) is a nilpotent matrix. Currently we test vanishing
// for powers under 1000 and fail with an error if none of these powers are zero.
//
// We follow the construction of p.72 of Lambek, Hofmann "Lectures on Rings and Modules", and
// we also follow their notation.

proc mfStrictifyIdempotent(matrix A, matrix U)
{
    int useSanityChecks = blowFlags("sanity_checks"); 
    
    // If U is already strictly idempotent, return it:
    if( U * U == U ){ return(U); }
    
    int m = ncols(A);
    
    matrix NN = U * U - U;
    matrix z[m][m];
    
    matrix powerOfN = NN;
    int vanishingOrder;
    
    for( vanishingOrder = 1; vanishingOrder < 1000; vanishingOrder++ )
    {
        //dbprint(printlevel, "[mfStrictifyIdempotent] Looking at order " + string(vanishingOrder));
        if( powerOfN == z ){ break; }
        
        powerOfN = NN * powerOfN;
    }
    
    if( powerOfN != z )
    {
        print("[mfStrictifyIdempotent] Exceeded bound looking for nilpotent, exiting.");
        return();
    }
    
    // vanishingOrder is the smallest integer positive such that N^vanishingOrder = 0.
    int i;
    matrix X[m][m];
    
    for(i = 1; i < vanishingOrder; i++)
    {
        X = X + (-1)^(i+1) * int(binomial(2 * i, i)) * power(NN,i);
    }
    
    X = 1/2 * X;
    
    matrix E = U + X * (unitmat(m) - 2 * U);
    
    // Sanity check
    if( useSanityChecks )
    {
        if( A * E != E * A )
        {
            print("[mfStrictifyIdempotent] E is not a morphism.");
            return();
        }
    
        if( E * E != E )
        {
            print("[mfStrictifyIdempotent] Failed to strictify idempotent.");
            return();
        }
    }
    
    return(E);
}

////////////////////////////////////////////////////////////////////
// computeGradingFromInjection
//
// Given a matrix F interpreted as a morphism of free modules, and an
// intvec g assigning degrees 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.

proc computeGradingFromInjection(matrix F, intvec g)
{   
    // Check that the sizes of g and F are compatible:
    if( size(g) != nrows(F) )
    {
        print("[computeGradingFromInjection] Sizes of map and grading vector not compatible, exiting.");
        return();
    }
    
    if( F == matrix(0,nrows(F),ncols(F)) )
    {
        print("[computeGradingFromInjection] The given map is zero, exiting.");
        return();
    }
    
    int i,j;
    intvec r;
    int ncolsF = ncols(F);
    int sizeg = size(g);
    
    for(i=1;i<=ncolsF;i++)
    {
        int currentDegree = 0;
        for(j=1;j<=sizeg;j++)
        {
            if( F[j,i] != 0 )
            {
               if( currentDegree != 0 && currentDegree != g[j] + 2 * deg(F[j,i]) )
               {
                   print("[computeGradingFromInjection] Encountered multiple nonzero degrees, exiting.");
                   return();
               }
               else
               {
                   currentDegree = g[j] + 2 * deg(F[j,i]);
               }
            }
        }
        
        r[i] = currentDegree;
    }
    
    return(r);
}

////////////////////////////////////////////////////////////////////
// mathomoglist
//
// Given a matrix, returns a list with the homogeneous components,
// the first entry in the list being the degree zero elements.

proc mathomoglist(matrix A)
{
    int i = 0;
    matrix Apart[nrows(A)][ncols(A)];
    list L = Apart;
    while(A != Apart)
    {
        matrix B[nrows(A)][ncols(A)] = jet(A,i) - jet(A,i-1); // components of degree i
        L[i+1] = B;
        Apart = Apart + B;
        kill B;
        i++;
    }
    
    return(L);
}

////////////////////////////////////////////////////////////////////
// matrixSplitIdempotent
//
// Given an idempotent matrix of polynomials, provides a splitting, i.e.
// we return a pair (F,G) such that F * G = 1 and G * F = E.

proc matrixSplitIdempotent(matrix E)
{
    int useSanityChecks = blowFlags("sanity_checks"); 
    
    if( useSanityChecks )
    {
        if( E == matrix(0,nrows(E),ncols(E)) )
        {
            print("[matrixSplitIdempotent] Passed a zero matrix, exiting.");
            return();
        }
        
        if( E * E != E )
        {
            print("[matrixSplitIdempotent] Passed a non-idempotent matrix, exiting.");
            return();
        }
    }
    
    list Elist = mathomoglist(E);
    
    // First we split the scalar matrix Elist[1] in the usual way (this is the degree zero part of E):
    module M = Elist[1];
    
    if( useSanityChecks )
    {
        // minbase can only work on a homogeneous module.
        if( !homog(M) )
        {
            print("[matrixSplitIdempotent] Image of idempotent is not homogeneous, exiting.");
            return();
        }
    }
                
    // Compute a minimal generating set for the image. Note that minbase
    // is currently bugged (Singular 3-1-1) and may not create a linearly
    // independent generating set on the run through.
    module im = minbase(M);
    kill M;
    
    // If the image is zero, our splitting is just a pair of zero matrices.
    // Otherwise, we use lift to compute the map factorising e(i) through its image.
    if( im[1] == [0] )
    {
        print("[matrixSplitIdempotent] Failed to split idempotent, unexpected zeros, exiting.");
        return();
    }
    
    module n = lift(im,Elist[1]);
    matrix g(0) = matrix(im,ncols(Elist[1]),size(im));
    matrix f(0) = matrix(n,size(im),ncols(Elist[1]));
    
    if( f(0) * g(0) != unitmat(nrows(f(0))) || g(0) * f(0) != Elist[1] )
    {
        print("[matrixSplitIdempotent] minbase and lift splitting failed, exiting.");
        return();
    }
    
    // The final F and G will have the same dimensions as f0 and g0.
    matrix F[nrows(f(0))][ncols(f(0))] = f(0);
    matrix G[nrows(g(0))][ncols(f(0))] = g(0);
    list prefixList = unitmat(nrows(g(0)));
    
    // The prefixes are for example (writing e(i) = Elist[i+1] for the degree i part of e):
    // prefixList[1] = 1
    // prefixList[2] = e(1)^2 + e(2)
    // prefixList[3] = e(3) + e(2)e(1) + e(1)e(2) + e(1)^3.
    
    // We keep modifying F and G until they work! (Or we hit a cutoff.)
    int i = 1;
    while( (F * G != unitmat(nrows(F)) || G * F != E) && i < 1000 )
    {
        if( size(Elist) <= i ){ Elist[i+1] = matrix(0,nrows(Elist[1]),ncols(Elist[1])); }
        
        matrix f(i) = f(0) * Elist[i+1];
        matrix prefix[nrows(Elist[1])][ncols(Elist[1])];
        
        int j;
        for(j=1;j<=i;j++)
        {
            prefix = prefix + Elist[i-j+2] * prefixList[j];
        } 
        prefixList = prefixList + list(prefix);
        
        matrix g(i) = prefix * g(0);
        kill prefix, j;
        
        // Update F and G in preparation for the while check:
        F = F + f(i);
        G = G + g(i);
        kill f(i), g(i);
        i++;
    }
    
    if( F * G != unitmat(nrows(F)) || G * F != E )
    {
        print("[matrixSplitIdempotent] Failed to split the idempotent in the end, exiting.");
        return();
    }
    
    // print("[matrixSplitIdempotent] Split idempotent matrix in " + string(i-1) + " steps.");
    
    return(list(F,G));
}

////////////////////////////////////////////////////////////////////
// mfSplitIdempotent
//
// Takes as input an odd block matrix A interpreted as the differential on a Z/2-graded free module X,
// and an even block matrix E interpreted as an endomorphism of X, satisfying E * E = E. We return
// a list B, F, G consisting of the differential B on a MF Y splitting E and morphisms F: X -> Y,
// G: Y -> X such that E = G * F and 1_Y = F * G.
//
// If we are given a third param we assume it is a grading vector for X and use it to compute
// a grading vector r on Y making F,G homogeneous of degree zero. Our return is then B,F,G,r.
//
// NOTE: There is a special case. If the potential is zero and the splitting is a Z/2-graded
// module with different ranks in each degree ("nonsquare matrices") then we return a list
//
//      "nonsquare",f(0),f(1),g(0),g(1)
// 
// where f(i), g(i) splits the idempotent in degree i (possibly f(i) * g(i) = 0 if the splitting
// in some degree is the zero module). If there is a grading, we return an additional pair r(0)
// r(1) giving the grading in each degree separately.

proc mfSplitIdempotent(matrix A, matrix E, list #)
{   
    int useSanityChecks = blowFlags("sanity_checks"); 
    int minbaseCutoff = blowFlags("minbase_cutoff");
    
    int m = ncols(A);
    
    int useGrading;
    if( size(#) > 0 )
    {
        intvec Agr = #[1];
        int useGrading = 1;
    }
    
    if( useSanityChecks )
    {
        if( A * E != E * A )
        {
            print("[mfSplitIdempotent] Passed a non-morphism, exiting.");
            return();
        }
    
        if( E * E != E )
        {
            print("[mfSplitIdempotent] Passed a non-idempotent endomorphism, exiting.");
            return();
        }
    }
    
    // print("[mfSplitIdempotent] Badness of worst coeff in E: " + string(size(string(worstMonomialCoeff(E)))));
        
    // Extract the odd and even part of the differential and morphism:
    list blocks = extractblockmat(A);
    list e_blocks = extractblockmat(E);
    
    matrix a(0) = blocks[3]; matrix a(1) = blocks[2];
    matrix e(0) = e_blocks[1]; matrix e(1) = e_blocks[4];
    
    // We produce morphisms f0, f1, g0, g1 where fi, gi is a splitting of ei, in the
    // sense that gi * fi = ei and fi * gi = 1.
    int hitZeros(0), hitZeros(1);
    
    int i;
    for(i=0; i<=1; i++)
    {
        if( e(i) != matrix(0,nrows(e(i)),ncols(e(i))) )
        {
            // First try the most awesome splitting method:
            def trySplit = matrixSplitIdempotent(e(i));
            
            if( typeof(trySplit) == "list" )
            {
                // We succeeded!
                matrix f(i) = trySplit[1];
                matrix g(i) = trySplit[2];
            }
            else
            {
                // The new way failed, try the old way:
                print("");
                print("[mfSplitIdempotent] New method failed [!######################!]");
                print("");
                
                // M is the image of the idempotent.
                module M = e(i);

                if( useSanityChecks )
                {
                    // minbase can only work on a homogeneous module
                    if( !homog(M) )
                    {
                        print("[mfSplitIdempotent] Image of idempotent is not homogeneous, exiting.");
                        return();
                    }
                }

                // Compute a minimal generating set for the image. Note that minbase
                // is currently bugged (Singular 3-1-1) and may not create a linearly
                // independent generating set on the run through. So we try reapplying
                // minbase a certain number of times (set by the flag minbaseCutoff).
                module im = minbase(M);
                kill M;

                // If the image is zero, our splitting is just a pair of zero matrices.
                // Otherwise, we use lift to compute the map factorising e(i) through its image.
                if( im[1] == [0] )
                {
                    print("[mfSplitIdempotent] Image of idempotent was zero, impossible, exiting.");
                    return();
                }

                module n = lift(im,e(i));
                matrix g(i) = matrix(im,ncols(e(i)),size(im));
                matrix f(i) = matrix(n,size(im),ncols(e(i)));

                // We are now guaranteed that g(i) * f(i) = e(i), and by construction
                // the image of g(i) is equal to the image of e(i) (the module M).
                // Minbase is supposed to guarantee that the columns of g(i) are a minimal
                // generating set for this image (suppose this is the case).
                //
                // The splitting of e(i) is free (it is projective, and projective modules
                // over polynomial rings are free) and is isomorphic to the image of e(i),
                // so the image of g(i) is a free module. But then any generating set of
                // this image must contain a basis, so a minimal generating set _is_ a
                // basis. That is, the columns of g(i) must be a basis for the image.
                // Equivalently, the map g(i) is a monomorphism. Since its image is
                // equal to the image of e(i), it follows that f(i) is an epimorphism.
                // Then in the usual way we deduce that f(i) * g(i) = 1.

                // So: failure of f(i) * g(i) = 1 (which does happen!) shows us that in
                // fact the columns of g(i) are not a minimal generating set. This is a
                // bug in minbase, and the workaround is as follows:

                int minbaseCount = 1;
                while( f(i) * g(i) != unitmat(nrows(f(i))) && minbaseCount < minbaseCutoff )
                {
                    // Try running minbase again:
                    dbprint(printlevel, "[mfSplitIdempotent] WARNING: Encountered minbase bug, running it again...");
                    im = minbase(im);
                    n = lift(im,e(i));
                    g(i) = matrix(im,ncols(e(i)),size(im));
                    f(i) = matrix(n,size(im),ncols(e(i)));

                    minbaseCount++;
                }
                kill n, im;
            }
        }
        else
        {
            hitZeros(i) = 1;
            matrix g(i) = matrix(0,ncols(e(i)),1);
            matrix f(i) = matrix(0,1,ncols(e(i)));
        }
    }
    
    if( useSanityChecks )
    {
        int i;
        for(i=0; i<=1; i++)
        {
            // Do g and f factorise the idempotent?
            if( g(i) * f(i) != e(i) )
            {
                print("[mfSplitIdempotent] Splitting was not successful, as g * f != e in degree " + string(i) + ", exiting.");
                return();
            }
            
            // We should also have g * f = 1. The only excuse is if we hit zeros.
            if( !hitZeros(i) && f(i) * g(i) != unitmat(nrows(f(i))) )
            {
                print("[mfSplitIdempotent] Splitting was not successful, as f * g != 1 in degree " + string(i) + ", exiting.");
                print( g(i) * f(i) );
                //print("[mfSplitIdempotent] The above is the matrix that should be 1, it has been exported as errMat.");
                return();
            }
            
            // If we did hit zeros, then f * g should be a 1x1 zero matrix.
            if( hitZeros(i) && f(i) * g(i) != matrix(0,1,1) )
            {
                print("[mfSplitIdempotent] Splitting was not successful, error in degree " + string(i));
                print("[mfSplitIdempotent] f: " + string(nrows(f(i))) + "x" + string(ncols(f(i))) + " and g: " + string(nrows(g(i))) + "x" + string(ncols(g(i))));
                return();
            }
        }
    }
    
    // Construct the gradings on each degree. If the splitting is zero
    // in a degree then the grading is just the intvec (0).
    if( useGrading )
    {
        intvec Bgr(0), Bgr(1);
        intvec Agr(0) = Agr[1..(m div 2)];
        intvec Agr(1) = Agr[(m div 2+1)..m];
        
        int i;
        for(i=0; i<=1; i++)
        {
            if( !hitZeros(i) )
            {
                Bgr(i) = computeGradingFromInjection(g(i), Agr(i));
            }
        }
    }
        
    // If we have nonsquare matrices we have to do something special:
    if( hitZeros(0) || hitZeros(1) || nrows(f(0)) != nrows(f(1)) )
    {
        // This only makes sense if the potential is zero, so warn and quit otherwise:
        if( useSanityChecks )
        {
            if( A * A != matrix(0,nrows(A),ncols(A)) )
            {
                print("[mfSplitIdempotent] Found nonsquare matrices for nonzero potential, exiting.");
                return();
            }
        }
        
        // The case of zero potential and non-square splitting is special, so our return
        // value is special. We put a string in the first position of the return list, so
        // that the receiver knows to interpret the rest of the output in the special manner.
        
        list rls = "nonsquare",f(0),f(1),g(0),g(1);
        kill f(0), f(1), g(0), g(1);
        
        if( useGrading )
        {   
            rls = rls + list(Bgr(0),Bgr(1));
        }
        
        return(rls);
    }

    // Put the f and g matrices back together:
    matrix zf1[nrows(f(0))][ncols(f(1))];
    matrix zf2[nrows(f(1))][ncols(f(0))];
    
    matrix zg1[nrows(g(0))][ncols(g(1))];
    matrix zg2[nrows(g(1))][ncols(g(0))];
    
    matrix F = blockmat_nsq(f(0), zf1, zf2, f(1));
    matrix G = blockmat_nsq(g(0), zg1, zg2, g(1));
    kill f(0), f(1), g(0), g(1), zf1, zf2, zg1, zg2;
    
    // The differential on the splitting is given by:
    matrix B = F * A * G;
    
    //print("[mfSplitIdempotent] Badness of worst coeff in final B: " + string(size(string(worstMonomialCoeff(B)))));
    //print("[mfSplitIdempotent] Worst coeff in final B: " + string(worstMonomialCoeff(B)));
    
    list retll = B,F,G;
    
    if( useGrading )
    {
        // We are dealing with graded MFs.
        intvec Bgr = Bgr(0), Bgr(1);
        retll = retll + list(Bgr);
        
        if( useSanityChecks && useGrading )
        {
            if( !isGradingValid(B, Bgr) )
            {
                print("[mfSplitIdempotent] Grading on splitting is not valid, exiting.");
                return();
            }
            
            // Check that F and G are degree zero:
            if( !checkMorphismDegree(A, Agr, B, Bgr, F, 0) || !checkMorphismDegree(B, Bgr, A, Agr, G, 0) )
            {
                print("[mfSplitIdempotent] One of the maps F,G is not degree zero, exiting.");
                return();
            }
        }
    }
    
    return(retll);
}

proc findConstant(matrix A)
{
    int i,j,q;
    int nc = ncols(A);
    int nr = nrows(A);
    int cutoff = 10;

    for(q=1;q<=cutoff;q++)
    {
        for(i=q;i<=nr;i=i+cutoff)
        {
            for(j=1;j<=nc;j++)
            {
                if( A[i,j] != 0 )
                {
                    if( deg(A[i,j]) == 0 )
                    {
                        intvec v = i,j;
                        return(list(1,v));
                    }
                }
            }
        }
    }
    
    return(list(0,"na"));
}

// Much slower than findConstant but will keep down the size of coefficients:
proc findConstantMin(matrix A)
{
    number small;
    number sizz;
    intvec v;
    
    int i,j;
    int nc = ncols(A);
    int nr = nrows(A);
    
    for(i=1;i<=nr;i++)
    {
        for(j=1;j<=nc;j++)
        {
            if( A[i,j] != 0 )
            {
                if( deg(A[i,j]) == 0 )
                {
                    number c = number(A[i,j]);

                    if( small == 0 || absValue(numerator(c)) * absValue(denominator(c)) < sizz )
                    {
                        v = i,j;
                        small = c;
                        sizz = absValue(numerator(small)) * absValue(denominator(small));
                    }
                    kill c;
                }
            }
        }
    }
    
    if( small == 0 )
    {
        return(list(0,"na"));
    }
    else
    {
        return(list(1,v));
    }
}

////////////////////////////////////////////////////////////////////
// contentMat
//

proc contentMat(matrix A)
{
    module M = A;
    int i;
    vector v;
    for(i=1;i<=size(M);i++)
    {
        v = v + content(M[i])*gen(i);
    }
    
    return(content(v));
}

////////////////////////////////////////////////////////////////////
// worstMonomialCoeff
//
// Returns the worst monomial coefficient in a matrix. Here "worst" means
// the reduced fraction a/b with a * b the largest.

proc worstMonomialCoeff(matrix A)
{
    number big;
    intvec w;
    
    int i,j, k;
    int nc = ncols(A);
    int nr = nrows(A);

    poly allVars = 1;
    for(i=1;i<=nvars(basering);i++)
    {
        allVars = allVars * var(i);
    }
    
    for(i=1;i<=nr;i++)
    {
        for(j=1;j<=nc;j++)
        {
            if( A[i,j] != 0 )
            {
                // Look at all coeffs of monomials:
                matrix co = coef(A[i,j], allVars);
                for(k=1;k<=ncols(co);k++)
                {
                    number c = number(co[2,k]);
                    if( big == 0 || absValue(numerator(c)) * absValue(denominator(c)) > absValue(numerator(big)) * absValue(denominator(big)) )
                    {
                        w = i,j;
                        big = c;
                    }
                    kill c;
                }
            }
        }
    }
    
    return(big);
}

////////////////////////////////////////////////////////////////////
// checkMorphismDegree
//
// Given the differentials A, B of two MFs and an even supermatrix F
// representing a morphism A -> B, together with grading vectors
// Agr and Bgr making A,B into graded MFs, we return 1 if F represents
// a degree d morphism and 0 otherwise.

proc checkMorphismDegree(matrix A, intvec Agr, matrix B, intvec Bgr, matrix F, int d)
{
    // Sanity check:
    if( F * A != B * F )
    {
        print("[isMorphismDegreeZero] Not passed valid morphism, exiting.");
        return(0);
    }
    
    matrix f0 = extractblockmat(F)[1];
    matrix f1 = extractblockmat(F)[4];
    int blockCols = ncols(F) div 2;
    int blockRows = nrows(F) div 2;
    
    int i, j;
    for(i = 1; i <= blockRows; i++) // row
    {
        for(j = 1; j <= blockCols; j++) // column
        {
            if( f0[i,j] != 0 )
            {
                if( d - Bgr[i] + Agr[j] != 2 * deg(f0[i,j]) )
                {
                    dbprint(printlevel, "[isMorphismDegreeZero] First test failed at (" + string(i) + "," + string(j) + ")");
                    return(0);
                }
            }
            
            if( f1[i,j] != 0 )
            {
                if( d - Bgr[blockRows+i] + Agr[blockCols+j] != 2 * deg(f1[i,j]) )
                {
                    dbprint(printlevel, "[isMorphismDegreeZero] Second test failed at (" + string(i) + "," + string(j) + ")");
                    return(0);
                }
            }
        }
    }
    
    return(1);
}

////////////////////////////////////////////////////////////////////
// isGradingValid
//
// Given the differential of a MF and a grading vector, checks whether
// the pair defines a graded MF in the sense explained in "Usage Guide" above.
// Returns 1 if the input is a valid graded MF, 0 otherwise. Note that we
// do not test MFs of 0, we just return 1.

proc isGradingValid(matrix A, intvec g)
{
    // Find the degree of our potential. Let R denote the basering, which
    // we view as a graded ring with the grading which assigns degree 2 to
    // each variable.
    
    int m = ncols(A) div 2;
    poly W = (A * A)[1,1];
    
    if( W == 0 ){ return(1); }
        
    // In R this W has degree 2 * (n+1) where n is given by
    int n = deg(W) - 1;
    
    matrix a0 = extractblockmat(A)[3];
    matrix a1 = extractblockmat(A)[2];
    
    // TODO: this doesn't actually check that the entries are homogeneous!
    int i, j;
    for(i = 1; i <= m; i++) // row
    {
        for(j = 1; j <= m; j++) // column
        {
            // First we check that the polynomial a0[i,j] is "valid". It is a map
            // R{a} -> R{b} where a = g[j] and b = g[m + i], and for A to be a valid
            // graded MF this map should be of degree n + 1 (note doubling of grading!).
            
            // Note that zeros are always OK
            if( a0[i,j] != 0 )
            {
                if( n+1 - g[m+i] + g[j] != 2 * deg(a0[i,j]) )
                {
                    dbprint(printlevel, "[isGradingValid] First test " + string(n+1-g[m+i]+g[j]) + "=" + string(2 * deg(a0[i,j])) + " failed for (" + string(i) + "," + string(j) + ")");
                    return(0);
                }
            }
            
            // Next check that a1[i,j] is valid:
            if( a1[i,j] != 0 )
            {
                if( n+1 - g[i] + g[m+j] != 2 * deg(a1[i,j]) )
                {
                    dbprint(printlevel, "[isGradingValid] Second test " + string(n+1-g[i]+g[m+j]) + "=" + string(2 * deg(a1[i,j])) + " failed for (" + string(i) + "," + string(j) + ")");
                    return(0);
                }
            }
        }
    }
    
    // If we are still here then A is valid.
    return(1);
}

////////////////////////////////////////////////////////////////////
// matrixreduce
//
// After some bad experiences with examples we cannot trust Singular's gaussian elimination
// routines, so here is a primitive version which is good enough for our purposes.
//
// The input is a square matrix M. The output is a matrix N and four invertible matrices P, Q
// formatted as a triple (N, P, Q, Pinv, Qinv) such that P M Q = N. If no entry of M is a unit
// then N = M and P = Q = id. Otherwise the top left entry of N is a 1 and the first column
// of N contains no other nonzero entries. Pinv = P^{-1} and Qinv = Q^{-1}.
//
// We also assume that v = i,j is an intvec with M[i,j] a nonzero constant.

proc matrixreduce(matrix M, intvec v)
{
    int mrank = ncols(M);
    int i = v[1];
    int j = v[2];

    if( M[i,j] == 0 || deg(M[i,j]) != 0 )
    {
        print("[matrixreduce] Passed coordinates of a zero, or nonconstant, matrix entry, exiting.");
        return(0);
    }
    
    matrix permr = permrow(unitmat(mrank),1,i);
    matrix permc = permcol(unitmat(mrank),1,j);
                
    // We found our constant term, permute it to (1,1):
    matrix P = 1/(M[i,j]) * permr;
    matrix Pinv = M[i,j] * permr;
    matrix Q = permc;
    matrix Qinv = permc;
                
    M = P * M * Q;

    // Go through the first column and transform away everything:
    matrix Pn = unitmat(mrank);
    matrix Pninv = Pn;
    
    int k;
    for(k=2;k<=mrank;k++)
    {
        if(M[k,1] != 0)
        {
            Pn[k,1] = -M[k,1];
            Pninv[k,1] = M[k,1];
        }
    }
    
    M = Pn * M;
    P = Pn * P;
    Pinv = Pinv * Pninv;
                
    list l = M,P,Q,Pinv,Qinv;
    kill P, Pinv, Q, Qinv, Pn, Pninv, permr, permc;
    return(l);
}

////////////////////////////////////////////////////////////////////
// MFtensor
//
// Takes as input a pair (A,B) of odd block 2x2 matrices which we interpret as a pair of differentials
// on Z/2-graded free modules X,Y resp. The return is the differential on the tensor product X x Y.

proc MFtensor(matrix A, matrix B)
{
    // There are two differentials A x 1 and B x 1 on the tensor product and we return their sum.
    int arank = ncols(A);
    int brank = ncols(B);
    
    return( ZZtensor( unitmat(arank), B ) + ZZtensor( A, unitmat(brank) ) );
}

proc MFtensorGradings(intvec g, intvec 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 graded free modules
    // with the gradings described by intvecs g, h resp. then the intvec describing
    // the shifts 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;
    
    intvec 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);
}

proc MFtensorBigradings(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 bigraded free modules
    // with the bigradings described by lists g, h resp. then the list describing
    // the shifts 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);
}

////////////////////////////////////////////////////////////////////
// MFdual
//
// Takes as input an odd block 2x2 matrix A which we interpret as a differential
// on a Z/2-graded free modules X. Returns the differential of the "dual" MF, namely
// if A = (0 a1 newline a0 0) then we return (0 transpose(a0) newline -transpose(a1) 0)

proc MFdual(matrix A)
{
    list blocks = extractblockmat(A);
    matrix z[nrows(A) div 2][ncols(A) div 2];
    matrix dual = blockmat_nsq( z, transpose(blocks[3]), -transpose(blocks[2]), z );
    return(dual);
}

proc MFdualGrading(intvec v)
{
    // We take the grading vector and multiply it by minus one.
    return(-v);
}

////////////////////////////////////////////////////////////////////
// ZZtensor 
//
// Takes as input a pair (A,B) of block matrices which we interpret as a pair of maps of Z/2-graded
// free modules A: X -> X', B: Y -> Y' according to the decompositions X = X0 + X1 and Y = Y0 + Y1,
// similarly for X', Y'. We assume that A, B are both homogeneous (i.e. either even or odd), but we
// allow X and X' (resp. Y and Y') to have different ranks, so A,B may not be square.
//
// Consider the tensor product Z = X x Y (resp. Z' = X' x Y') of Z/2-graded modules and the map
// A x B: Z -> Z' which sends a homogeneous tensor f x g to (-1)^{|f||B|} A(f) x B(g). We return
// the matrix of this map with respect to the decomposition Z0 = (X0 x Y0) + (X1 + Y1) and
// Z1 = (X0 + Y1) x (X1 + Y0) of Z (resp. of Z'). The basis of X0 x Y0 etc. is ordered in the
// usual way for a Kronecker product of matrices.

proc ZZtensor(matrix A, matrix B)
{
    list la = extractblockmat(A);
    list lb = extractblockmat(B);
    
    matrix a11, a12, a21, a22, b11, b12, b21, b22;
    
    a11 = la[1]; a12 = la[2]; a21 = la[3]; a22 = la[4];
    b11 = lb[1]; b12 = lb[2]; b21 = lb[3]; b22 = lb[4];
    
    int degb = 0; // Default B to degree 0, i.e. even
    
    matrix zero[nrows(B) div 2][ncols(B) div 2];
    matrix bigzero[nrows(B)][ncols(B)];
    
    if( B != bigzero && b11 == zero && b22 == zero ){ degb = 1; }
    
    // Define the four corners of Z, then use blockmat:
    matrix z11 = blockmat_nsq( tensor(a11,b11), (-1)^(degb) * tensor(a12,b12), tensor(a21,b21), (-1)^(degb) * tensor(a22,b22) );
    matrix z12 = blockmat_nsq( tensor(a11,b12), (-1)^(degb) * tensor(a12,b11), tensor(a21,b22), (-1)^(degb) * tensor(a22,b21) );
    matrix z21 = blockmat_nsq( tensor(a11,b21), (-1)^(degb) * tensor(a12,b22), tensor(a21,b11), (-1)^(degb) * tensor(a22,b12) );
    matrix z22 = blockmat_nsq( tensor(a11,b22), (-1)^(degb) * tensor(a12,b21), tensor(a21,b12), (-1)^(degb) * tensor(a22,b11) );
        
    matrix Z = blockmat_nsq(z11,z12,z21,z22);
    
    return(Z);
}

////////////////////////////////////////////////////////////////////
// extractblockmat takes a block matrix (A B // C D) and returns a list {A, B, C, D}.
// We assume that the blocks each have half the width of the matrix, and half the height
// of the matrix (so we assume both are even, but we allow nonsquare blocks).

proc extractblockmat(matrix block)
{
    int nr = nrows(block) div 2;
    int nc = ncols(block) div 2;
    
    matrix a[nr][nc];
    matrix b[nr][nc];
    matrix c[nr][nc];
    matrix d[nr][nc];
    
    int i,j; // i is the row and j the column
    for(i=1;i <= nr; i++)
    {
        for(j=1;j <= nc; j++)
        {
            a[i,j] = block[i,j];
            b[i,j] = block[i,j+nc];
            c[i,j] = block[i+nr,j];
            d[i,j] = block[i+nr,j+nc];
        }
    }
    
    list l = a,b,c,d;
    kill a,b,c,d,nr,nc,i,j;
    return(l);
}

////////////////////////////////////////////////////////////////////
// fuseDefects
//
// Returns the reduced tensor product of Y with X, as a pair (B, e) 
// consisting of the differential B of a matrix factorisation of the total
// potential W and an idempotent endomorphism e whose splitting is (YxX)[n],
// where n is the number of internal variables.

proc fuseDefects(matrix Y, matrix X, poly W)
{
    int prl = printlevel;

    // First compute the internal algebra. That is, find powers of the internal
    // variables which can be written as linear combinations of the partials. 

    list t = varTrans(W);
    matrix C = t[2];
    intvec q = t[1];
    
    // The ideal whose quotient gives the algebra J is generated by the internal
    // variables to the powers given by q.
    int i;
    ideal J;
    int numX = numXVars();
    
    for( i = 1; i <= numX; i++ )
    {
        J[i] = x(i)^(q[i]);
    }

    // Compute the ordinary tensor product:
    matrix D = MFtensor(Y,X);
    
    // Inflate or blow up the differential on this tensor product, that is,
    // compute the differential of Y x X x I.
    dbprint(prl, "[fuseDefects] Blowing up differential on tensor product.");
    matrix Dblow = mablow(D,J);
    
    // Reduce this inflated differential:
    dbprint(prl, "[fuseDefects] Reducing blown up differential.");
    list l = mfReduce(Dblow);
    matrix RD = l[1];
    matrix F = l[2]; // A homotopy equivalence F: Dblow -> RD
    matrix G = l[3]; // A homotopy equivalence G: RD -> Dblow

    // Compute the matrix dQ_X:
    dbprint(prl, "[fuseDefects] Computing dQ.");
    matrix dq = dQ(X, q, C);
    matrix iddQ = ZZtensor( unitmat(ncols(Y)), dq ); // "iddQ" is "id x dQ_X"
    
    // Compute the matrix deltaQ:
    dbprint(prl, "[fuseDefects] Computing deltaQ.");
    matrix delQ = deltaQ(D, q);

    // And finally the idempotent e on Dblow:
    int sign = (numX * (numX-1)) div 2;
    matrix e = (-1)^(sign) * iddQ * delQ;
    
    // The induced idempotent on the reduction is denoted ep:
    matrix ep = F * e * G;

    // 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("[fuseDefects] The reduced idempotent is not a morphism, exiting.");
        return();
    }
    
    list l = RD, ep;
    
    return(l);
}


proc supertrace(matrix m)
{
    int r = ncols(m) div 2;
    matrix zero[r][r];
    matrix E[r][r]; E = E + 1; // The unit matrix

    matrix B = blockmat(E,zero,zero,-1 * E);

    return( trace( B * m ) );
}

proc mfReduce(matrix A, list #)
{
    int useSanityChecks = blowFlags("sanity_checks");
    int factorOutContent = blowFlags("factor_out_content");
    
    dbprint(printlevel, "[mfReduce] Reducing matrix of size " + string(ncols(A)));
        
    system("--ticks-per-sec",1000);
    int timeElapsed = timer;
    
    int useGrading;
    if( size(#) > 0 )
    {
        intvec Agr = #[1];
        int useGrading = 1;
    }
    
    if( factorOutContent )
    {
        // Pull out the content of the matrix, to reduce to the size of coeffs being multiplied:
        number Acont = contentMat(A);
        //print("[mfReduce] Content of A: " + string(Acont));
        A = A * 1/Acont;
    }

    // Extract the odd and even part of the differential:
    list blocks = extractblockmat(A);
    matrix b0 = blocks[3]; matrix b1 = blocks[2];
    kill blocks;

    // Reduce b0. This returns a list b0r,P,Q,n where P,Q are invertible matrices,
    // b0r is a matrix with no nonzero constant entries, and P * b0 * Q = dsum(unitmat(n),b0r).
    // (Unless n = 0 in which case P * b0 * Q = b0r.) 
    // The induced differential in degree one is dsum(W * unitmat(n),b1r) = Q^{-1} * b1 * P^{-1}.
    
    list l = zoomreduce(b0);
    // n = l[2];
    // P = l[3];
    // Q = l[4];
    // Pinv = l[5];
    // Qinv = l[6];

    b0 = l[1];
    b1 = submat(l[6] * b1 * l[5],(l[2]+1)..nrows(b1),(l[2]+1)..ncols(b1));
    
    matrix f0 = submat(l[6],(l[2]+1)..nrows(l[6]),1..ncols(l[6]));
    matrix f1 = submat(l[3],(l[2]+1)..nrows(l[3]),1..ncols(l[3]));
    matrix g0 = submat(l[4],1..nrows(l[4]),(l[2]+1)..ncols(l[4]));
    matrix g1 = submat(l[5],1..nrows(l[5]),(l[2]+1)..ncols(l[5]));

    if( useGrading )
    {
        // The grading on P * b0 * Q, Q^{-1} * b1 * P is induced from Q and Pinv.
        intvec Agr0 = Agr[1..(ncols(A) div 2)];
        intvec Agr1 = Agr[(ncols(A) div 2 + 1)..ncols(A)];
        
        Agr0 = computeGradingFromInjection(l[4], Agr0);
        Agr1 = computeGradingFromInjection(l[5], Agr1);
        
        // Next in both differentials we split off the first n columns and n rows,
        // so do the same to the grading:
        Agr0 = Agr0[(l[2]+1)..size(Agr0)];
        Agr1 = Agr1[(l[2]+1)..size(Agr1)];
    }
    
    kill l;
        
    // Now repeat the awesome feat on b1. This time b1 is replaced by P * b1 * Q
    // and b0 is replaced by Qinv * b0 * Pinv (with a nxn trivial piece removed).
    list l = zoomreduce(b1);

    b1 = l[1];
    b0 = submat(l[6] * b0 * l[5],(l[2]+1)..nrows(b0),(l[2]+1)..ncols(b0));
    
    f0 = submat(l[3] * f0,(l[2]+1)..nrows(l[3]),1..ncols(f0));
    f1 = submat(l[6] * f1,(l[2]+1)..nrows(l[6]),1..ncols(f1));
    g0 = submat(g0 * l[5],1..nrows(g0),(l[2]+1)..ncols(l[5]));
    g1 = submat(g1 * l[4],1..nrows(g1),(l[2]+1)..ncols(l[4]));
    
    if( useGrading )
    {
        // The grading on Qinv * b0 * Pinv, P * b1 * Q is induced from Pinv and Q.
        Agr0 = computeGradingFromInjection(l[5], Agr0);
        Agr1 = computeGradingFromInjection(l[4], Agr1);
        
        // Next in both differentials we split off the first n columns and n rows,
        // so do the same to the grading:
        Agr0 = Agr0[(l[2]+1)..size(Agr0)];
        Agr1 = Agr1[(l[2]+1)..size(Agr1)];
    }
    
    kill l;

    // The maps F = (f0,f1) and G = (g0,g1) now give homotopy equivalences between
    // the original MF and the reduced MF (b0,b1). 
    matrix zb[nrows(b0)][ncols(b0)];
    matrix zf[nrows(f0)][ncols(f0)];
    matrix zg[nrows(g0)][ncols(g0)];

    matrix B = blockmat(zb, b1, b0, zb);
    matrix F = blockmat_nsq(f0, zf, zf, f1);
    matrix G = blockmat_nsq(g0, zg, zg, g1);
    if( useGrading ){ intvec Bgr = Agr0, Agr1; }
    kill zb,zf,zg,b1,b0,f0,f1,g0,g1;
    
    // Partial check that F * G is homotopic to 1. Note that the
    // check we use won't work to check G * F simeq 1 because we
    // need the MF in question to be reduced.
    if( useSanityChecks )
    {
        // Check that F and G are morphisms
        if( F * A != B * F || G * B != A * G )
        {
            print("[mfReduce] Either F or G is not a morphism, exiting.");
            return();
        }
        
        if( useGrading )
        {
            if( !isGradingValid(B, Bgr) )
            {
                print("[mfReduce] Produced a non-valid grading on reduction, exiting.");
                return();
            }
        }
        
        if( 0 ) // debug additional checks
        {
            matrix shouldBeZero = F * G - unitmat(nrows(F));
            matrix Bc = B;
            
            int j;
            for(j=1; j<=nvars(basering); j++)
            {
                shouldBeZero = subst(shouldBeZero,var(j),0);
                Bc = subst(Bc,var(j),0);
            }
        
            if( shouldBeZero != matrix(0,nrows(shouldBeZero),ncols(shouldBeZero)) )
            {
                print("[mfReduce] F * G is not homotopic to the identity, exiting.");
                return();
            }
            
            // Check that B is actually reduced
            if( Bc != matrix(0,nrows(Bc),ncols(Bc)) )
            {
                print("[mfReduce] Error, the return matrix is not reduced, exiting.");
                return();
            }
        }
    }

    //dbprint(printlevel, "[mfReduce]          worst output coeff: " + string(worstMonomialCoeff(B)));
    dbprint(printlevel, "[mfReduce]          result is size " + string(ncols(B)) + ".");
    dbprint(printlevel, "[mfReduce]          elapsed time " + string(timer-timeElapsed) + "ms.");
    
    if( factorOutContent )
    {
        B = B * Acont;
    }
    
    list ret = B, F, G;
    kill B, F, G;
    
    if( useGrading )
    {
        ret = ret + list(Bgr);
    }
    
    return(ret);
}

////////////////////////////////////////////////////////////////////
// isHomotopyIdempotent
//
// Partial check whether the given matrix is an idempotent up to homotopy.

proc isHomotopyIdempotent(matrix E)
{
    int ycount = numYVars();
    int i;
    for(i=1;i<=ycount;i++)
    {
        E = subst( E, y(i), 0 );
    }
    
    if( E * E != E )
    {
        return(0);
    }
    
    return(1);
}

////////////////////////////////////////////////////////////////////
// mfPushforward
//
// We are given a matrix A interpreted as the differential on a matrix factorisation X
// of a potential W in our ring, together with a ring variable "intvar" and integer N
// such that intvar^N acts null-homotopically on X, and a matrix H of degree one such that
// A * H + H * A = intvar^N * unitmat, i.e. the entries of H give a homotopy between
// intvar^N and the zero map. We assume that, as a polynomial, W does not depend on the
// variable intvar.
//
// Let S denote the ring obtained from our basering by deleting intvar. We can restrict
// scalars on X and obtain an infinite rank MF of W over S, the so-called "pushforward".
// The first return value of this routine is the differential of a finite rank MF Y of W over S
// homotopy equivalent to this pushforward.
//
// Let Xblow denote the MF of W over S obtained by "inflation". Then the pushforward is a
// direct summand of Xblow[1] and the second and third return value are morphisms F: Xblow[1] -> Y
// and G: Y -> Xblow[1] such that F * G == unitmat and G * F is, up to homotopy, the same
// idempotent on Xblow[1] determined by its direct summand the pushforward.
//
// If we are passed an intvec grading for A as the last param, we respect this grading
// throughout and return as a fourth param a grading vector on Y, such that Y is equivalent
// as a graded MF to the pushforward.
//
// We assume that our ring is R = 0,(some variables, we don't care),dp.
//
// NOTE: Before calling any of the routines of blow.lib we rename intvar to x(1) and all the other
// variables to y's, so that our active ring is of the form 0,(x(1),y(1),...,y(??)),dp and
// S is 0,(y(1),...,y(??)),dp. Naturally we convert back to the original ring before returning
// our results.
//
// Note: if the fourth param is a _list_ of matrices then we assume that the first is the
// H described above, and the rest are homotopies on A. We then return an additional list consisting
// of the homotopies induced on Y by these homotopies in the list.
//
// NOTE: If the potential is zero and splitting an idempotent produces nonsquare matrices, then 
// the return format is different, because the Z/2-graded module may be zero in one degree: this
// means we can't write things as block matrices. In this case we return a list 
//
//  "nonsquare", f0,f1,g0,g1 [possibly two grading vectors r1,r2]
//
// where f0,g0 split the idempotent on Xblow[1] in degree zero and resp. for f1,g1. Hence the
// differential on Y is implicit (obtained from fi,gi and the differential on Xblow[1]). The
// grading vectors ri give the degrees on the even and odd part of Y resp.

proc mfPushforward(matrix A, poly intvar, int N, Hlist, list #)
{
    int useSanityChecks = blowFlags("sanity_checks"); 
    dbprint(printlevel, "[mfPushforward] Pushing forward matrix of size " + string(ncols(A)) + " with N = " + string(N));
    //dbprint(printlevel, "[mfPushforward] Worst input coeff: " + string(worstMonomialCoeff(A)));
    
    system("--ticks-per-sec",1000);
    int timeElapsed = timer;
    
    def RRR = basering;
    
    // Check for grading:
    int useGrading, returnAblow;
    if( size(#) > 0 )
    {
        if( typeof(#[1]) == "intvec" )
        {
            intvec Agr = #[1];
            useGrading = 1;
        
            if( useSanityChecks )
            {
                if( !isGradingValid(A,Agr) )
                {
                    print("[mfPushforward] Given an invalid grading, exiting.");
                    return();
                }
            }
        }
        
        if( typeof(#[1]) == "string" )
        {
            // At the end of our return list, add Ablow
            if( #[1] == "return_ablow" )
            {
                returnAblow = 1;
            }
        }
    }
    
    // Check for a list of homotopies to track:
    int homotopyList;
    if( typeof(Hlist) == "list" )
    {
        matrix H = Hlist[1];
        Hlist = delete(Hlist,1);    
        
        if( size(Hlist) != 0 )
        {
            int homotopyList = 1;
        }
    }
    else
    {
        matrix H = Hlist;
    }
    
    if( useSanityChecks )
    {
        if( A * H + H * A != intvar^N * unitmat(ncols(A)) )
        {
            print("[mfPushforward] Not passed valid homotopy, exiting.");
            return();
        }
    }
    
    // Convert everything into a more suitable ring. We go through the original
    // list of ring variables, rename everything apart from intvar to y(i)'s (in
    // ascending order) and rename intvar to x(1).
    list rlist = ringlist(RRR);
    list varlist = rlist[2];
    
    // Create the new list of variables:
    list newvar;
    int i;
    int ycount;
    for(i=1;i<=size(varlist);i++)
    {
        if( string(intvar) == varlist[i] )
        {
            newvar = newvar + list("x(1)");
        }
        else
        {
            ycount++;
            string s = "y(" + string(ycount) + ")";
            newvar = newvar + list(s);
            kill s;
        }
    }
    kill varlist;
    
    list newringList;
    newringList[1] = rlist[1];
    newringList[2] = newvar;
    newringList[3] = rlist[3];
    newringList[4] = rlist[4];
    kill rlist;
    
    // Fix the variable weighting
    intvec kk = (1..size(newvar));
    for(i=1; i<=size(newvar); i++)
    {
        kk[i] = 1;
    }
    newringList[3][1][2] = kk;
    kill kk;

    int useParam;
    if( defined(minpolyblow) && defined(r) ) // deal with minpoly
    { 
        newringList[1][4][1] = 0;
        useParam = 1;
    }
    
    // Create our new ring:
    def nR = ring(newringList);
    setring nR;
    
    // Now complete the correct definition of nR by specifying the right minpoly: 
    if( useParam )
    {
        // Note that fetch uses the position of the ring variables, not the name, so 
        // that fetch(var(1)) is always var(1).
        poly minpolyblow = fetch(RRR,minpolyblow); 
        poly zz = subst(minpolyblow,var(1),r);
        number nu = leadcoef(zz);
        minpoly = nu;
        poly minpolyblow = fetch(RRR,minpolyblow);
        export(minpolyblow);
        // NOTE: The double definition of minpolyblow is because Singular nukes
        // all variables in our ring when we set minpoly
    }    
    
    // If we call fetch then the object of RRR is converted to the new renamed variables
    matrix A = fetch(RRR,A);
    matrix H = fetch(RRR,H);
    if( homotopyList )
    {
        list Hlist = fetch(RRR,Hlist);
    }
    
    // So now H * A + A * H = x(1)^N * unitmat and the potential W depends only on the y's.
    // Compute the inflation:	
    ideal J = x(1)^N;
    matrix Ablow = mablow(A,J);
    if( useGrading ){ intvec Ablowgr = mablowGrading( Agr, J, N ); }

    if( homotopyList )
    {
        for(i=1;i<=size(Hlist);i++)
        {
            Hlist[i] = mablow(Hlist[i],J);
        }
    }
    
    // Compute the idempotent on the inflation:
    ideal J = x(1)^N; 
    matrix e = mablow( H, J ) * mablow_delta(A, 1, N, J);
    	
    // Reduce the inflated matrix factorisation (of W over S = QQ[y(1),...,y(??)]):
    if( !useGrading )
    {
        // for mfReduceAlt use list l = mfReduceAlt(Ablow,e).
        list l = mfReduce(Ablow);
    }
    else
    {
        // For mfReduceAlt use list l = mfReduceAlt(Ablow, e, Ablowgr).
        list l = mfReduce(Ablow, Ablowgr);
    }
    matrix RT = l[1]; // Differential on reduced MF
    matrix F = l[2]; // A homotopy equivalence Ablow -> RT
    matrix G = l[3]; // A homotopy equivalence RT -> Ablow
    
    // For mfReduceAlt use matrix ep = l[useGrading + 4].
    matrix ep = F * e * G;
    
    if( useGrading )
    {
        intvec RTgr = l[4];
    }
    
    kill l;

    // Move the homotopies over to the reduction RT:
    if( homotopyList )
    {
        for(i=1;i<=size(Hlist);i++)
        {
            Hlist[i] = F * Hlist[i] * G;
        }
    }
    
    if( useSanityChecks )
    {
        if( ep * RT != RT * ep )
        {
            print("[mfPushforward] The reduced idempotent (up to homotopy) ep is not an endomorphism of RT, exiting.");
            return();
        }
    }
    
    // Check that ep is idempotent up to homotopy, and strictify it:
    matrix epstrict;
    if( ep * ep != ep )
    {
        matrix epck = ep;
        for(i=1;i<=ycount;i++)
        {
            epck = subst( epck, y(i), 0 );
        }

        if( epck * epck != epck )
        {
            // NOTE: If this test fails to pass, it is usually due to a sign problem...
            if( epck * epck == -1 * epck )
            {
                ep = -ep;
            }
        }
        kill epck;

        epstrict = mfStrictifyIdempotent(RT, ep);
    }
    else
    {
        epstrict = ep;
    }
    
    if( useSanityChecks )
    {
        if( epstrict * RT != RT * epstrict || epstrict * epstrict != epstrict )
        {
            print("[mfPushforward] The strictified idempotent is not an endomorphism of RT, or not an idempotent on the nose, exiting.");
            return();
        }
        
        if( useGrading )
        {
            if( !checkMorphismDegree( RT, RTgr, RT, RTgr, epstrict, 0 ) )
            {
                print("[mfPushforward] Strictified idempotent is not degree zero, exiting.");
                return();
            }
        }
    }
    
    list l;
    if( !useGrading )
    {
        l = mfSplitIdempotent(RT, epstrict);
    }
    else
    {
        l = mfSplitIdempotent(RT, epstrict, RTgr);
    }

    // This can only happen if the potential is zero and mfSplitIdempotent found non-square
    // matrices. In this case our return formatting is different, see the introduction.
    if( typeof(l[1]) == "string" )
    {
        // In this the list l has the format l = "nonsquare",f0,f1,g0,g1 where fi,gi are the
        // splittings in the respective degrees, and l may have as a final entry a grading
        // vector for the even and odd parts (separately).
        
        // IMPORTANT: We need to "suspend" this answer in order to get the splitting
        // which is equal to Y. 
        
        setring RRR;
        list l = imap(nR,l);
        list retlist = l;
        retlist[2] = l[3]; // suspend by replacing f0 by f1
        retlist[3] = l[2]; // and f1 by f0
        retlist[4] = l[5];
        retlist[5] = l[4];
        
        if( useGrading )
        {
            retlist[6] = l[7]; // also suspend the gradings
            retlist[7] = l[6];
        }
        kill l;
        
        dbprint(printlevel, "[mfPushforward] elapsed time " + string(timer-timeElapsed) + "ms.");
            
        return(retlist);
    }
    
    // mfSplitIdempotent returns not just the splitting of the idempotent, but the
    // actual maps splitting the idempotent. In the following notation, this means
    // that s * t == ep and t * s == unitmat. Here t = l[2] and s = l[3].
    
    matrix final = l[1];
    if( useGrading ){ intvec finalgr = l[4]; }
    
    matrix t = l[2]; // t: RT --> final
    matrix s = l[3]; // s: final --> RT
    
    // We have now split the original idempotent e on Ablow, and the splittings are:
    matrix mapAblowToFinal = t * F;
    matrix mapFinalToAblow = G * s;
    
    // Each of the homotopies in our list is a homotopy on RT, let us replace
    // these by homotopies on final.
    if( homotopyList )
    {
        for(i=1;i<=size(Hlist);i++)
        {
            Hlist[i] = t * Hlist[i] * s;
        }
    }
    
    kill s,t,F,G,RT,e,ep,epstrict;
    if( !returnAblow ){ kill Ablow; }
    if( useGrading ){ kill RTgr; }

    // Back to the original ring and return. Note that we use fetch and not imap
    // because we want to replace variables in nR by the variables in RRR in the same
    // _position_ in ringlist.
    setring RRR;
    matrix final = fetch(nR,final);
    matrix mapAblowToFinal = fetch(nR,mapAblowToFinal);
    matrix mapFinalToAblow = fetch(nR,mapFinalToAblow);
    if( returnAblow ){ matrix Ablow = fetch(nR, Ablow); }
    if( homotopyList )
    {
        list Hlist = fetch(nR,Hlist);
    }
    
    // Note that we really want a splitting of e[1] on Ablow[1], so let us suspend:
    final = mfSuspend( final );
    mapAblowToFinal = mfSuspendMorph( mapAblowToFinal );
    mapFinalToAblow = mfSuspendMorph( mapFinalToAblow );
    if( useGrading ){ finalgr = mfSuspendGrading( finalgr ); }
    if( returnAblow ){ Ablow = mfSuspend( Ablow ); }
    
    if( homotopyList )
    {
        for(i=1;i<=size(Hlist);i++)
        {
            Hlist[i] = mfSuspend( Hlist[i] ); // homotopies are suspended just like differentials
        }
    }
    
    if( !returnAblow )
    {
        // Normal behaviour:
        list retlist = list(final,mapAblowToFinal,mapFinalToAblow);
    }
    else
    {
        // Special behaviour:
        list retlist = list(final, mapAblowToFinal,mapFinalToAblow, Ablow);
    }
    
    // Note that grading (if it exists) goes before homotopy list (if it exists!).
    if( useGrading )
    {
        retlist = retlist + list(finalgr);
    }

    if( homotopyList )
    {
        retlist = retlist + list(Hlist);
        kill Hlist;
    }
    
    dbprint(printlevel, "[mfPushforward] result is size " + string(ncols(final)) + ".");
    dbprint(printlevel, "[mfPushforward] elapsed time " + string(timer-timeElapsed) + "ms.");
    kill final, mapAblowToFinal, mapFinalToAblow;
        
    return(retlist);
}

////////////////////////////////////////////////////////////////////
// mfPushforwardInductive
//
// We are given a matrix A interpreted as the differential on a matrix factorisation X
// of a potential W in our ring, together with a list of ring variables "intvars" and list NN
// such that intvar[i]^NN[i] acts null-homotopically on X, and a matrix H[i] of degree one such that
// A * H[i] + H[i] * A = intvar[i]^NN[i] * unitmat, i.e. the entries of H[i] give a homotopy between
// intvar[i]^NN[i] and the zero map. We assume that, as a polynomial, W does not depend on the
// variables intvar.
//
// Let S denote the ring obtained from our basering by deleting the intvars. We can restrict
// scalars on X and obtain an infinite rank MF of W over S, the so-called "pushforward".
// The first return value of this routine is the differential of a finite rank MF Y of W over S
// homotopy equivalent to this pushforward.
//
// The second return value is the grading vector for A, if we are using gradings.
//
// The third (or second, if there is no grading) return value is a list R (of the same length
// as intvars) whose entries are tuples
//
//       (F,G,u,NN,D[,r])
//
// which records the F and G matrices returned by mfPushforward at each stage of the computation,
// along with the internal variable u in question and the exponent NN it is associated with. The D
// is the differential of the pushforward, and r is its optional grading vector. Note that if in
// the splitting we encounter non-square matrices the format is slightly different; see the code.

proc mfPushforwardInductive(matrix A, list intvars, list NN, list H, list #)
{
    int useSanityChecks = blowFlags("sanity_checks"); // Turn off to run faster
    
    system("--ticks-per-sec",1000);
    int timeElapsed = timer;
    
    def RRR = basering;
    list R;
    
    // Check for grading:
    int useGrading;
    if( size(#) > 0 )
    {
        intvec Agr = #[1];
        useGrading = 1;
    }
    
    // If we are not passed any internal variables, we perform no pushforward:
    if( size(intvars) == 0 )
    {
        list retlist = A;
        if( useGrading )
        {
            retlist = retlist + list(Agr);
        }
        return(retlist);
    }
    
    int i;
    if( useSanityChecks )
    {
        for(i=1;i<=size(intvars);i++)
        {
            if( A * H[i] + H[i] * A != intvars[i]^(NN[i]) * unitmat(ncols(A)) )
            {
                print("[mfPushforward] Not passed valid homotopies, exiting.");
                return();
            }
        }
    }
    
    // Inductively we contract off the action of each of the external variables, using
    // the given homotopies. During the process we switch rings, killing off the 
    // external variables that we have used up.
    
    int count = 1;
    def nR(1) = basering;
    int useParam;
    
    while( size(intvars) > 0 )
    {
        dbprint(printlevel, "[mfPushforwardInductive] In step " + string(count) + " with ring vars " + varstr(nR(count)));
        
        // Do the reduction with respect to the first variable and homotopy, and
        // induce homotopies on this reduction for the action of the rest of the variables.
        if( !useGrading )
        {
            list l = mfPushforward(A, intvars[1], NN[1], H);
        }
        else
        {
            list l = mfPushforward(A, intvars[1], NN[1], H, Agr);
        }
        
        if( typeof(l[1]) == "string" )
        {   
            // mfSplitIdempotent found, during mfPushforward, that it had to split
            // an idempotent on a complex. Our return value in this situation is
            // "nonsquare",[r1,r2],R' where r1,r2 are possible gradings. We are given
            // "nonsquare",f0,f1,g0,g1,[r1,r2] by mfPushforward.
            poly ivar = intvars[1];
            setring RRR;
            
            if( count > 1 )
            {
                list l = imap(nR(count),l);
                poly ivar = imap(nR(count),ivar);
            }
            
            list ll = "nonsquare";
            
            if( useGrading )
            {
                ll = ll + list(l[6],l[7]); // Add on r1,r2
            }
            
            // Modify the current map list R by adding on the end the list
            //
            //    ("nonsquare",f0,f1,g0,g1,intvars[1],NN[1][,r1,r2])
            //
            // NOTE: We can't stick these together into block matrices (like the other
            // entries in R) because we might be rank zero in one component
            
            list tempp = "nonsquare",l[2],l[3],l[4],l[5],ivar,NN[1];
            if( useGrading ){ tempp = tempp + list(l[6],l[7]); }
            R = R + list(tempp);
            kill tempp;
            
            ll = ll + list(R);
                        
            dbprint(printlevel, "[mfPushforwardInductive] total elapsed time " + string(timer-timeElapsed) + "ms.");
            
            return(ll);
        }
        
        // Legend
        // final = l[1] is the MF with intvars[1] excluded which is h.e. to the pushforward
        // l[2] is a map Ablow[1] --> final
        // l[3] is a map final --> Ablow[1] (where [1] means suspension)
        // l[4] is the grading (if we are using gradings) or our modified list of homotopies
        // l[5] is our homotopies if we are using the grading, and undefined otherwise
    
        // Return to the original ring and add to R the tuple (l[2],l[3],intvar[1],NN[1],l[1][,l[4]]). Note
        // we have to return to the original ring because l[2],l[3] may depend on
        // external variables which we are going to delete in a moment.
        poly ivar = intvars[1];
        
        setring RRR;
        
            // The first time through nR(1) = RRR so there is no need for imap.
            if( count > 1 )
            {
                list l = imap(nR(count),l);
                poly ivar = imap(nR(count),ivar);
            }

            // Add the tuple (F,G,u,NN,D[,r])
            list tempp = l[2], l[3], ivar, NN[1], l[1];
            if( useGrading ){ tempp = tempp + list(l[4]); }
            R = R + list(tempp);
            kill tempp;
            
        setring nR(count);
        
        // Remove the variable intvars[1] from our active ring:
        list rlist = ringlist(nR(count));
        list varlist = rlist[2];
    
        // Create the new list of variables:
        int i;
        for(i=1;i<=size(varlist);i++)
        {
            if( varlist[i] == string(intvars[1]) )
            {
                varlist = delete(varlist,i);
                break;
            }
        }
    
        list newringList;
        newringList[1] = rlist[1];
        newringList[2] = varlist;
        newringList[3] = rlist[3];
        newringList[4] = rlist[4];
        // Fix the variable weighting
        intvec kk = (1..size(varlist));
        for(i=1; i<=size(varlist); i++)
        {
            kk[i] = 1;
        }
        newringList[3][1][2] = kk;
        
        if( defined(minpolyblow) && defined(r) ) // deal with minpoly
        { 
            newringList[1][4][1] = 0;
            useParam = 1;
        } 
        
        // Prepare for the next step of the induction:
        A = l[1];
        if( useGrading ){ Agr = l[4]; }
        
        intvars = delete(intvars,1);
        NN = delete(NN,1);
        
        if( size(intvars) == 0 )
        {
            break;
        }

        H = l[4 + useGrading];
    
        // Create our new ring:
        def nR(count+1) = ring(newringList);
        setring nR(count+1);
    
        // Now complete the correct definition of nR by specifying the right minpoly: 
        if( useParam )
        {
            // Note that fetch uses the position of the ring variables, not the name, so 
            // that fetch(var(1)) is always var(1).
            poly minpolyblow = fetch(nR(count),minpolyblow); 
            poly zz = subst(minpolyblow,var(1),r);
            number nu = leadcoef(zz);
            minpoly = nu;
			poly minpolyblow = fetch(nR(count),minpolyblow);
	        export(minpolyblow);
    	    // NOTE: The double definition of minpolyblow is because Singular nukes
        	// all variables in our ring when we set minpoly
        }
        
        matrix A = imap(nR(count),A);
        list intvars = imap(nR(count),intvars);
        list H = imap(nR(count),H);
        
        count++;
    }
    
    // Return to the original ring:
    setring RRR;
    matrix Anew = imap(nR(count),A);
    
    list retlist = Anew;
    
    if( useGrading )
    {
        retlist = retlist + list(Agr);
    }
    
    // Add on the maps:
    retlist = retlist + list(R);
    dbprint(printlevel, "[mfPushforwardInductive] total elapsed time " + string(timer-timeElapsed) + "ms.");
    
    return(retlist);
}

proc displaySplittingRecord(list R)
{
    print("[displaySplittingRecord] Given a splitting record of length " + string(size(R)));
    
    int i;
    for(i=1;i<=size(R);i++)
    {
        if( size(R[i]) == 4 )
        {
            string out = "[displaySplittingRecord]    Variable: ";
            out = out + string(R[i][3]) + " exponent: " + string(R[i][4]);
            out = out + " F: " + string(nrows(R[i][1])) + "x" + string(ncols(R[i][1])) + " G: ";
            out = out + string(nrows(R[i][2])) + "x" + string(ncols(R[i][2]));
            print(out);
        }
        else
        {
            string out = "[displaySplittingRecord]    Nonsquare, variable: ";
            out = out + string(R[i][5]) + " exponent: ";
            out = out + string(R[i][6]) + " f0: " + string(nrows(R[i][1])) + "x" + string(ncols(R[i][1]));
            out = out + " f1: " + string(nrows(R[i][2])) + "x" + string(ncols(R[i][2]));
            out = out + " g0: " + string(nrows(R[i][3])) + "x" + string(ncols(R[i][3]));
            out = out + " g1: " + string(nrows(R[i][4])) + "x" + string(ncols(R[i][4]));
            print(out);
        }// f0,f1,g0,g1,intvars[1],N[1]
    }
}

////////////////////////////////////////////////////////////////////
// isMorphismNull
//
// Given two MFs QQQ, PPP (in block format) and a morphism F: QQQ --> PPP
// returns 1 if the map is null-homotopic and zero otherwise.

proc isMorphismNull(matrix QQQ, matrix PPP, matrix F)
{
    if( PPP * F != F * QQQ )
    {
        print("[isMorphismNull] Given matrix is not a morphism, exiting.");
        return();
    }
    
    vector B = matrix2vector(F);
    
    matrix f = blocksq(QQQ)[3];
    matrix g = blocksq(QQQ)[2];
    matrix fp = blocksq(PPP)[3];
    matrix gp = blocksq(PPP)[2];
    
    module M=syz(dminus(f,g,fp,gp));
    matrix Mker=M;
    
    // Test that B (= F as a vector) actually belongs to the kernel of dminus:
    if( reduce(B, std(Mker)) != [0] )
    {
        print("[isMorphismNull] Strange error, exiting.");
        return();
    }
    
    // Test whether B is in the image of dplus
    module im = module(dplus(f,g,fp,gp));
    
    if( reduce(B, std(im)) == [0] )
    {
        return(1);
    }
    else
    {
        return(0);
    }
}

////////////////////////////////////////////////////////////////////
// chernChar
//
// Given a matrix factorisation, returns the Chern character as a
// reduced element modulo the Jacobi ideal (so two MFs have the same
// chern char, as an element of the Jacobi algebra, iff the return
// values of this routine are the same). By definition
//
//      ch(A) = (-1)^{n choose 2} str(diff(A,x(1)),...,diff(A,x(n))
// 
// where n is the total number of variables.

proc chernChar(matrix A)
{
    // First determine the potential:
    poly W = (A*A)[1,1];
    
    // Now the Jacobi algebra
    ideal jaci = jacob(W);
    ideal J = std(jaci);
    int numVars = nvars(basering);

    int i;
    matrix X = unitmat(ncols(A));
    
    for(i=1; i<=numVars;i++)
    {
        X = X * diff(A,var(i));
    }
    
    int p = (numVars*(numVars-1)) div 2;
    poly ch = (-1)^(p) * supertrace(X);
    ch = reduce( ch, J );
    return(ch);
}

////////////////////////////////////////////////////////////////////
// Old routines for computing cohomology
//
// Let Q = (0 g newline f 0) and Q' = (0 g' newline f' 0) be matrix factorisations
// of a common potential.
// 
// We define the operator dminus on even elements of Hom(Q,Q') (the Z2
// graded Hom, before taking cohomology) and dplus on odd elements, by
//
// dminus(a) = Q' a - a Q    and   dplus(a) = Q' a + a Q.
//
// i.e. the differential on Hom(Q,Q'). To compute cohomology we want to
// find the kernel and images of these maps, but first we convert them
// into matrices. We do this by writing dminus (resp. dplus) as 2rr' x 2rr'
// matrices which act on maps a represented as column vectors (in the natural
// way, see the preprint). This is done in the routines dminus and dplus
// below. Then we can take Ker and Im and, using matrix2vector and vector2matrix,
// convert between the representations.

////////////////////////////////////////////////////////////////////
// rdsum
//
// Returns the matrix A^{oplus r}

proc rdsum(matrix A, int r)
{
    list B;
    B[1]=A;
    for (int i=2; i<=r; i++)
    {
        B[i]=dsum(B[i-1],A);
    }
    return(B[r]);
}

proc dminus(matrix f,g,f',g')
{
    int r=nrows(f);
    int r'=nrows(f');
    matrix P1=concat(tensor(f',unitmat(r)), rdsum(-transpose(f),r'));
    matrix P2=concat(rdsum(transpose(-g),r'),tensor(g',unitmat(r)));
    matrix P[2*r*r'][2*r*r']=P2,P1;
    return(P);
}

proc dplus(matrix f,g,f',g')
{ 
    int r=nrows(f);
    int r'=nrows(f');
    matrix P1=concat(rdsum(transpose(f),r'),tensor(g',unitmat(r)));
    matrix P2=concat(tensor(f',unitmat(r)),rdsum(transpose(g),r'));
    matrix P[2*r*r'][2*r*r']=P1,P2;
    return(P);
}

proc vector2matrix(list V, int r, int r')
{
    int i,j,k;
    list M;
    matrix ma[2*r'][2*r];
    for(i=1; i<=size(V); i++){M = M + list(ma);}
    for(i=1; i<=size(V); i++)
    {
        for(j=1; j<=r; j++)
        {
            for(k=1; k<=r'; k++)
            {
                M[i][k,j] = V[i][(k-1)*r + j, 1];
                M[i][k+r',j+r] = V[i][(r'-1+k)*r + j, 1];
            }
        }
    }

    return(M);
}

proc vector2matrixFerm(list V, int r, int r')
{
    int i,j,k;
    list M;
    matrix ma[2*r'][2*r];
    for(i=1; i<=size(V); i++){M = M + list(ma);}
    for(i=1; i<=size(V); i++)
    {
        for(j=1; j<=r; j++)
        {
            for(k=1; k<=r'; k++)
            {
                M[i][k,j+r] = V[i][(k-1)*r + j, 1];
                M[i][k+r',j] = V[i][(r'-1+k)*r + j, 1];
            }
        }
    }
    
    return(M);
}

proc matrix2vector(matrix M)
{
    int r1 = ncols(M) div 2;
    int r2 = nrows(M) div 2;
    int i,ii,j,k;
    vector v;
    
    // bosonic part:
    for(ii=1; ii<= r2; ii++)
    {
        for(j=1; j<= r1; j++)
        {
            v = v + gen((ii-1)*r1 + j)*M[ii,j];
        }
    }

    for(ii=1; ii<= r2; ii++)
    {
        for(j=1; j<= r1; j++)
        {
            v = v + gen(r1*r2 + (ii-1)*r1 + j)*M[ii+r2,j+r1];
        }
    }

    // fermionic part:
    for(i=1; i<= r2; i++)
    {
        for(j=1; j<= r1; j++)
        {
            v = v + gen(2*r1*r2 + (i-1)*r1 + j)*M[i,j+r1];
        }
    }

    for(i=1; i<= r2; i++)
    {
        for(j=1; j<= r1; j++)
        {
            v = v + gen(2*r1*r2 + r1*r2 + (i-1)*r1 + j)*M[i+r2,j];
        }
    }

    return(v);
}

////////////////////////////////////////////////////////////////////
// MFco
//
// We are given the differentials (0 g newline f 0) and (0 g' newline f' 0) of two
// matrix factorisations X,Y and we return a basis over the basefield for
// the cohomology H of the two-periodic complex Hom(X,Y), as a list
// MM with MM[1] a basis of H^0 and MM[2] a basis of H^1.

proc MFco(matrix f,g,f',g')
{
    module M=syz(dminus(f,g,f',g'));
    matrix Mker=M;
    module Q=modulo(M,module(dplus(f,g,f',g')));
    matrix base=kbase(std(Q));
    matrix m = Mker*base;

    module Mf=syz(dplus(f,g,f',g'));
    matrix Mkerf=Mf;
    module Qf=modulo(Mf,module(dminus(f,g,f',g')));
    matrix basef=kbase(std(Qf));
    matrix mf = Mkerf*basef;

    int i;
    int r = ncols(f);
    int r' = ncols(f');
    matrix ma[2*r'][2*r];
    list L;
    for(i=1; i<=ncols(m); i++){L = L + list(ma);}
    intvec v1 = 1..nrows(m);
    intvec v2;

    list Lf;
    for(i=1; i<=ncols(mf); i++){Lf = Lf + list(ma);}
    intvec vf1 = 1..nrows(mf);
    intvec vf2;

    for(i=1; i<=ncols(m); i++)
    {
        intvec v2 = i;
        L[i] = submat(m,v1,v2);
    }

    for(i=1; i<=ncols(mf); i++)
    {
        intvec vf2 = i;
        Lf[i] = submat(mf,vf1,vf2);
    }

    list MMb = vector2matrix(L, r, r');
    list MMf = vector2matrixFerm(Lf, r, r');

    matrix null[2*r'][2*r];
    list leer;
    if(MMb[1] == null)
        {MMb = leer;}
    if(MMf[1] == null)
        {MMf = leer;}

    list MM;
    MM[1] = MMb;
    MM[2] = MMf;

    return(MM);
}

proc MFcofull(matrix QQQ,PPP)
{
    matrix f = blocksq(QQQ)[3];
    matrix g = blocksq(QQQ)[2];
    matrix f' = blocksq(PPP)[3];
    matrix g' = blocksq(PPP)[2];

    list MM = MFco(f,g,f',g');
    return(MM);
}


// Old routines:

proc zoomreduce(matrix A)
{
    int useSmallCoeffs = blowFlags("small_coeffs");
    int useSanityChecks = blowFlags("sanity_checks");

    // n is the current size of the top left identity block, begins at zero
    int n = 0;
    int nr = nrows(A);
    int nc = ncols(A);
    int i;
    int ncolsOrig = ncols(A);
    if( useSanityChecks ){ matrix Aorig = A; }

    matrix P = unitmat(nr);
    matrix Q = unitmat(nc);
    matrix Pinv = unitmat(nr);
    matrix Qinv = unitmat(nc);

    // Always the final matrix B will be P * A * Q
    while( 1 )
    {
        if( ncolsOrig > 250 && (n mod 5) == 0 && n > 0 )
        {
            dbprint(printlevel, "[zoomreduce]    reduced to size " + string(2*(ncolsOrig - n)));
        }
        
        // Find a unit in the given matrix, transform it to the top left, and eliminate the first
        // column and first row
        
        // r[2][1] is the row and r[2][2] the column our found constant entry of A
        if( useSmallCoeffs )
        {
            list r = findConstantMin(A);
        }
        else
        {
            list r = findConstant(A);
        }

        if( r[1] != 1 ){ break; }
        
        if( r[2][1] != 1 )
        {
            A = permrow(A,1,r[2][1]);
            P = permrow(P,1+n,r[2][1]+n);
            Pinv = permcol(Pinv,1+n,r[2][1]+n);
        }

        if( r[2][2] != 1 )
        {
            A = permcol(A,1,r[2][2]);
            Q = permcol(Q,1+n,r[2][2]+n);
            Qinv = permrow(Qinv,1+n,r[2][2]+n);
        }
        kill r;

        // Note you have to do these types of things on P,Q first! Because A will change!
        if( A[1,1] != 1 )
        {
            P = multrow(P,1+n,1/A[1,1]);    
            Pinv = multcol(Pinv,1+n,A[1,1]);
            A = multrow(A,1,1/A[1,1]);
        }

        for(i=2;i<=nc;i++)
        {
            if( A[1,i] != 0 )
            {
                Q = addcol(Q,1+n,-A[1,i],i+n);
                Qinv = addrow(Qinv,i+n,A[1,i],1+n);
                A = addcol(A,1,-A[1,i],i);
            }
        }

        for(i=2;i<=nr;i++)
        {
            if( A[i,1] != 0 )
            {
                P = addrow(P,1+n,-A[i,1],i+n);
                Pinv = addcol(Pinv,i+n,A[i,1],1+n);
                A = addrow(A,1,-A[i,1],i);
            }
        }

        if( nr == 1 || nc == 1 ){ break; }
        
        A = submat(A,2..nr,2..nc);
        nc = nc - 1;
        nr = nr - 1;
        n = n + 1;
    }

    if( useSanityChecks ) //debug very slow
    {
        matrix Afin = A;
        if( n > 0 ){ Afin = dsum(unitmat(n),Afin); }
        
        if( Afin != P * Aorig * Q || P * Pinv != unitmat(ncols(P)) || Q * Qinv != unitmat(ncols(Q)) )
        {
            print("[zoomreduce] Error in reduction, exiting.");
            return();
        }
        
        kill Afin;
    }
    
    // The final reduced matrix is the direct sum of the identity matrix and our final A,
    // but we just return the interesting bit A. So P * original A * Q = dsum(unitmat(n),A)
    // and both P and Q are invertible
    return(list(A,n,P,Q,Pinv,Qinv));
}

////////////////////////////////////////////////////////////////////
// mfReduce
//
// Takes as input an odd block matrix A interpreted as the differential on a Z/2-graded free module X.
// The output is a list B, F, G. Here B is the differential of a matrix factorisation Y homotopy
// equivalent to X which has no constant entries (the "reduced" form of the original matrix
// factorisation). Note the size of the input and output matrices may differ.
//
// The second and third outputs are matrices F,G representing morphisms F: X -> Y and G: Y -> X resp.
// which are mutually inverse in the homotopy category.
//
// If we are optionally passed a grading intvec, then we return a list B, F, G, Bgr where Bgr
// is the grading vector on Y such that F, G are morphisms of degree zero.

proc mfReduceAlt(matrix A, matrix e, list #)
{
    int useSanityChecks = blowFlags("sanity_checks");
    dbprint(printlevel, "[mfReduceAlt] Reducing matrix of size " + string(ncols(A)));
    
    system("--ticks-per-sec",1000);
    int timeElapsed = timer;
    
    int useGrading;
    if( size(#) > 0 )
    {
        intvec Agr = #[1];
        intvec Agr0 = Agr[1..(ncols(A) div 2)];
        intvec Agr1 = Agr[(ncols(A) div 2 + 1)..ncols(A)];
        int useGrading = 1;
    }
    
    // Extract the odd and even part of the differential
    list blocks = extractblockmat(A);
    matrix b0 = blocks[3]; matrix b1 = blocks[2];
    kill blocks;
    
    matrix f0 = unitmat(ncols(b0));
    matrix f1 = unitmat(ncols(b1));
    matrix g0 = unitmat(ncols(b0));
    matrix g1 = unitmat(ncols(b0));
    matrix e0 = extractblockmat(e)[1];
    matrix e1 = extractblockmat(e)[4];
    
    // Throughout this whole process b0, b1 is a MF of W,
    // and (f0,f1) is a map A --> B, and (g0,g1) is a map B --> A,
    // and indeed these maps are homotopy inverse to one another.
    // Step by step we reduce the size of the matrices b0,b1 by peeling
    // off contractible summands. At each stage F is a split epi and G
    // is a split mono in each degree.
    
    // Find a unit in b1, transform it to the top left, and eliminate the first
    // column and first row using mfReduce_splitworker
    
    // r[2][1] is the row and r[2][2] the column our found constant entry of A
    // r[1] is 1 if we found a constant, and zero otherwise
    list r = findConstantMin(b1);
    int madeStrictIdempotent = 0;
    
    while( r[1] == 1 && !madeStrictIdempotent )
    {
        // First permute this constant to row 1 and column 1 of b1, and do
        // the corresponding transformation on b0 and the F,G maps.
        if( r[2][1] != 1 )
        {
            b1 = permrow(b1,1,r[2][1]);
            b0 = permcol(b0,1,r[2][1]);
            f0 = permrow(f0,1,r[2][1]);
            g0 = permcol(g0,1,r[2][1]);
            e0 = permrow(e0,1,r[2][1]);
            e0 = permcol(e0,1,r[2][1]);
        }

        if( r[2][2] != 1 )
        {
            b1 = permcol(b1,1,r[2][2]);
            b0 = permrow(b0,1,r[2][2]);
            f1 = permrow(f1,1,r[2][2]);
            g1 = permcol(g1,1,r[2][2]);
            e1 = permrow(e1,1,r[2][2]);
            e1 = permcol(e1,1,r[2][2]);
        }
        
        // Now the pair b0, b1 is suitable for feeding into mfReduce_splitworker
        list spL = mfReduce_splitworker(b0, b1);
        
        // The return list from this is p0,p1,u0,u1,v0,v1
        b0 = spL[1];
        b1 = spL[2];
        f0 = spL[3] * f0;
        f1 = spL[4] * f1;
        g0 = g0 * spL[5];
        g1 = g1 * spL[6];
        e0 = spL[3] * e0 * spL[5];
        e1 = spL[4] * e1 * spL[6];
        
        if( e0 * e0 == e0 && e1 * e1 == e1 )
        {
            madeStrictIdempotent = 1;
        }
        //else
        //{
        //    if( (isHomotopyIdempotent(e0) || isHomotopyIdempotent(-e0) ) &&
        //        (isHomotopyIdempotent(e1) || isHomotopyIdempotent(-e1) ) )
        //    {
        //        dbprint(printlevel, "[mfReduceAlt] Found early idempotent up to homotopy...");
        //        madeStrictIdempotent = 1;
        //    }
        //}
        
        kill spL;
        
        // This completes the step, now look for another constant
        if( !madeStrictIdempotent ){ r = findConstantMin(b1); }
    }
    
    // Now perform the same feat on b0
    if( !madeStrictIdempotent )
    {
        list r = findConstantMin(b0);

        while( r[1] == 1 && !madeStrictIdempotent )
        {
            if( r[2][1] != 1 )
            {
                b0 = permrow(b0,1,r[2][1]);
                b1 = permcol(b1,1,r[2][1]);
                f1 = permrow(f1,1,r[2][1]);
                g1 = permcol(g1,1,r[2][1]);
                e1 = permrow(e1,1,r[2][1]);
                e1 = permcol(e1,1,r[2][1]);
            }

            if( r[2][2] != 1 )
            {
                b0 = permcol(b0,1,r[2][2]);
                b1 = permrow(b1,1,r[2][2]);
                f0 = permrow(f0,1,r[2][2]);
                g0 = permcol(g0,1,r[2][2]);
                e0 = permrow(e0,1,r[2][2]);
                e0 = permcol(e0,1,r[2][2]);
            }

            // Now the pair b1, b0 is suitable for feeding into mfReduce_splitworker
            list spL = mfReduce_splitworker(b1, b0);

            // The return list from this is p1,p0,u1,u0,v1,v0
            b0 = spL[2];
            b1 = spL[1];
            f0 = spL[4] * f0;
            f1 = spL[3] * f1;
            g0 = g0 * spL[6];
            g1 = g1 * spL[5];
            e0 = spL[4] * e0 * spL[6];
            e1 = spL[3] * e1 * spL[5];

            if( e0 * e0 == e0 && e1 * e1 == e1 )
            {
                madeStrictIdempotent = 1;
            }
            //else
            //{
            //    if( (isHomotopyIdempotent(e0) || isHomotopyIdempotent(-e0) ) &&
            //        (isHomotopyIdempotent(e1) || isHomotopyIdempotent(-e1) ) )
            //    {
            //        dbprint(printlevel, "[mfReduceAlt] Found early idempotent up to homotopy...");
            //        madeStrictIdempotent = 1;
            //    }
            //}

            kill spL;

            // This completes the step, now look for another constant
            if( !madeStrictIdempotent ){ r = findConstantMin(b0); }
        }
    }

    // The final F and G are split epi and split mono resp., so we can use G
    // to determine the grading on B from the grading on A.
    if( useGrading )
    {
        intvec Bgr0 = computeGradingFromInjection(g0, Agr0);
        intvec Bgr1 = computeGradingFromInjection(g1, Agr1);
    }
    
    // The maps F = (f0,f1) and G = (g0,g1) now give homotopy equivalences between
    // the original MF and the reduced MF (b0,b1)
    matrix zb[nrows(b0)][ncols(b0)];
    matrix zf[nrows(f0)][ncols(f0)];
    matrix zg[nrows(g0)][ncols(g0)];

    matrix B = blockmat(zb, b1, b0, zb);
    matrix F = blockmat_nsq(f0, zf, zf, f1);
    matrix G = blockmat_nsq(g0, zg, zg, g1);
    matrix E = F * e * G;
    
    if( useGrading ){ intvec Bgr = Bgr0, Bgr1; }
    kill zb,zf,zg,b1,b0,f0,f1,g0,g1;
    
    dbprint(printlevel, "[mfReduceAlt]          result is size " + string(ncols(B)) + ".");
    dbprint(printlevel, "[mfReduceAlt]          elapsed time " + string(timer-timeElapsed) + "ms.");
        
    list ret = B, F, G;
    kill B, F, G;
    
    if( useGrading )
    {
        ret = ret + list(Bgr);
    }
    
    ret = ret + list(E);
    
    return(ret);
}

proc mfReduce_splitworker(matrix q0, matrix q1)
{
    // We are given two matrices q0, q1 of the same dimensions
    if( ncols(q0) < 2 || nrows(q0) < 2 || ncols(q1) != ncols(q0) || nrows(q1) != nrows(q0) || ncols(q0) != nrows(q0) )
    {
        print("[mfReduce_splitworker] Incorrect data, exiting.");
        return();
    }
    
    // We break q0 as (a,b newline c,d) and q1 as (e, f newline g,h ) and we assume
    // that q0 * q1 = q1 * q0 = W for some potential W, and that e is a
    // scalar.
    
    int n = ncols(q0);
    
    if( q1[1,1] == 0 || deg(q1[1,1]) != 0 )
    {
        print("[mfReduce_splitworker] e failed to be a unit, exiting.");
        return();
    }
    
    number e = number(q1[1,1]);
    
    matrix d = submat(q0,2..nrows(q0),2..ncols(q0));
    matrix h = submat(q1,2..nrows(q1),2..ncols(q1));
    matrix f = submat(q1,1,2..ncols(q1));
    matrix g = submat(q1,2..nrows(q1),1);
    
    // Then the pair (d, h - 1/egf) is a matrix factorisation of W and
    // it is homotopy equivalent to (q0,q1) via the maps
    //
    //    U = (u0,u1): (q0,q1) --> (d, h - 1/e gf)
    //    V = (v0,v1): (d, h - 1/egf) --> (q0,q1)
    // 
    // defined by u0 = (-1/eg, 1), u1 = (0 1), v0 = (0 1)^t, v1 = (-1/e f, 1)^t
    
    matrix p0 = d;
    matrix p1 = h - 1/e * (g * f);
    matrix u0 = concat(-1/e * g, unitmat(n-1));
    matrix u1 = concat(matrix(0,n-1,1),unitmat(n-1));
    matrix v0 = transpose(u1);
    matrix v1 = transpose(concat(-1/e * transpose(f),unitmat(n-1)));

    list retlist = p0,p1,u0,u1,v0,v1;
    kill e, d, h, f, g, p0, p1, u0, u1, v0, v1;
    return(retlist);
}

// Functions added since first paper post


/////////////////////////////////////////////////
// variableDeflationList

proc variableDeflationList(intvec v)
{
    int numX = numXVars();
    int i;
    list L;
    
    for( i=1; i<=numX; i++ )
    {
        L = L + list(varDeflate(v, i));
    }
    
    return(L);
}

/////////////////////////////////////////////////
// varDeflate
//
// Given an intvec (v1,...,vn) of the same length as
// the number of ring variables, returns the matrix
// which gives the action of multiplication by var(z)
// on the ring as a module over the subring generated
// by the ring variables to their respective powers, i.e.
// var(1)^v1, ..., var(n)^vn.
//
// BUT in the matrix of this action, which is a polynomial
// in these powers of the variables, we replace var(i)^vi
// by var(i) (i.e. we have restricted to the subring).

proc varDeflate(intvec v, int z)
{
    ideal J;
    int i;
    for( i=1; i<=size(v); i++)
    {
        if( v[i] <= 0 )
        {
            print("[varDeflate] Passed bad integers, exiting.");
        }
        
        J = J + var(i)^(v[i]);
    }

    ideal basis = kbase(std(J));
    int jdim = size(basis);
    
    list basisList;
    for( i=1; i<=jdim; i++)
    {
        basisList[i] = basis[i];
    }
    
    kill basis, J;
    
    // Now basis will be the monomials of var(i)-degree < v[i] in 
    // the ordering dictated by the ring. We find the matrix of the
    // action of var(z) on RRR as a module over the poly ring with powers of variables.
    
    list MM;
    
    for(i=1; i<=jdim; i++)
    {
        matrix M[jdim][1];
        poly f = var(z) * basisList[i];
        int s = indexinlist(basisList, f);

        if( s != 0 )
        {
            // f is still a basis element
            M[s,1] = 1;
        }
        else
        {
            // f is not a basis element, so we must have been 
            // var(z)^(v[z]-1) times a basis element.
            poly g = f/(var(z)^(v[z]));
        
            int t = indexinlist(basisList, g);    
            M[t,1] = var(z); // if this were var(z)^(v[z]) you would get the honest matrix
        }
        
        MM = MM + list(M);
    }
    
    matrix U = MM[1];
    for( i=2; i<=size(MM); i++)
    {
        U = concat(U, MM[i]);
    }

    return(U);
}

/////////////////////////////////////////////////
// indexinlist
//
// Given a list and anything p, returns 0 if p is
// not in l and the index into the list of p otherwise.

proc indexinlist(list l, p)
{
    int i;
    for(int i=1; i<=size(l); i++)
    {
        if(p==l[i])
        {
            return(i);
        }
    }
    return(0);
}

proc madefl(matrix M, intvec v)
{   
    // TODO export this globally
    list variableActions = variableDeflationList(v);
    int n = ncols(variableActions[1]);
    
    // Define L to be an appropriately indexed list of blown-up matrices:
    int i1,j1,i2,j2,i;
    
    int ncolsM = ncols(M);
    int nrowsM = nrows(M);
    
    list e,L;
    for(i=1; i<=nrowsM; i++)
    {
        L[i] = e;
    }
    for(i1=1; i1<=nrowsM; i1++) // row
    {
        for(j1=1; j1<=ncolsM; j1++) // column
        {
            matrix PM = poly2matrix(M[i1,j1], variableActions);
            L[i1][j1] = PM;
            kill PM;
        }
    }

    matrix A[nrowsM*n][ncolsM*n];

    for(i1=1; i1<=nrowsM; i1++)
    {
        for(j1=1; j1<=n; j1++)
        {
            for(i2=1; i2<=ncolsM; i2++)
            {
                for(j2=1; j2<=n; j2++)
                {
                    A[(i1-1)*n + j1, (i2-1)*n + j2] = L[i1][i2][j1,j2];
                }
            }
        }
    }

    return(A);
}