linkhom.lib

version="1.0";
category="Miscellaneous";
info="
LIBRARY:   linkhom.lib  Computation of link homology using web compilation
AUTHOR:    Nils Carqueville, Daniel Murfet
KEYWORDS:  matrix factorisation, link homology

PROCEDURES:
";

LIB "linalg.lib";
LIB "matrix.lib";
LIB "control.lib"; // for rightKernel
LIB "ring.lib";
LIB "blow.lib";
LIB "mfweb.lib";

////////////////////////////////////////////////////////////////////
// USAGE GUIDE
//
// We expect the ambient ring to have the same number of variables as edges
// in the link under consideration. In principle params should be OK, but
// this has not been checked.
//
// The references to [KR] are to Khovanov and Rozansky, "Matrix factorizations and
// link homology", arXiv:math/0401268v2.

////////////////////////////////////////////////////////////////////
// NOTES ON GRADING
//
// In [KR, p.47] all variables are assigned degree two. The conventions for graded
// matrix factorisations are explained in [KR, p.44]

/////////////////////////////////////////
// Link datatype
//
// Suppose we are given a link with m crossings, on which we fix a particular ordering.
// We represent this as a list (X,Y) consisting of an intvec X and a list Y, such that
//
// - X is an intvec of size m with X[i] = 1 if the i-th crossing of the link is positive
// and X[i] = -1 if the crossing is negative. We follow the positive/negative conventions of [KR].
// - each entry of Y is an intvec of size 4 such that Y[i] labels the edges incident at
//  the i-th crossing of the link, oriented in the usual clockwise fashion (see below).

////////////////////////////////////////////////////////////////////
// linkLibrary
//
// Returns various standard links, referenced by name.
//
// The ordering on edge labellings follows
//
//       a        b
//         \   /
//          \/
//         /\
//       /   \
//     d       c
//
// i.e. this crossing would be represented by a,b,c,d.

proc linkLibrary(string s)
{
    list retLink;
    
    // 3 crossing knots ////////////////////////////////////////
    
    if( s == "trefoil" || s == "3_1" )
    {
        // This is the right handed trefoil
        intvec edgevec(1) = 2,1,3,4;
        intvec edgevec(2) = 4,3,5,6;
        intvec edgevec(3) = 6,5,1,2;
        intvec vertexParity = -1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3));    
    }
    
    // 4 crossing knots ////////////////////////////////////////
    
    if( s == "figure-eight" || s == "4_1" )
    {
        intvec edgevec(1) = 1,7,8,5;
        intvec edgevec(2) = 6,8,7,2;
        intvec edgevec(3) = 3,5,6,4;
        intvec edgevec(4) = 4,2,1,3;
        intvec vertexParity = 1,1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4));
    }
    
    // 5 crossing knots ////////////////////////////////////////
    
    if( s == "cinquefoil" || s == "5_1" )
    {
        intvec edgevec(1) = 1,2,10,9;
        intvec edgevec(2) = 3,4,2,1;
        intvec edgevec(3) = 6,5,4,3;
        intvec edgevec(4) = 8,7,5,6;
        intvec edgevec(5) = 9,10,7,8;
        intvec vertexParity = 1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5));
    }
    
    if( s == "three-twist" || s == "5_2" )
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 5,6,2,1;
        intvec edgevec(3) = 4,10,9,5;
        intvec edgevec(4) = 7,9,10,8;
        intvec edgevec(5) = 8,3,6,7;
        intvec vertexParity = 1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5));
    }
        
    // 6 crossing knots ////////////////////////////////////////
    
    if( s == "stevedore" || s == "6_1" )
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 8,6,1,7;
        intvec edgevec(3) = 7,9,10,8;
        intvec edgevec(4) = 11,10,9,12;
        intvec edgevec(5) = 12,4,5,11;
        intvec edgevec(6) = 5,3,2,6;
        intvec vertexParity = -1,1,1,1,1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }
    
    if( s == "stevedore_alternate" || s == "6_1_alternate" )
    {
        intvec edgevec(1) = 12,1,2,3;
        intvec edgevec(2) = 4,2,1,11;
        intvec edgevec(3) = 6,3,4,5;
        intvec edgevec(4) = 5,7,8,6;
        intvec edgevec(5) = 10,8,7,9;
        intvec edgevec(6) = 9,11,12,10;
        intvec vertexParity = 1,1,-1,-1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }
    
    if( s == "6_1_bad" ) // Old 6_1, preserved for inspection
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 6,8,7,1; // This appears to be an error, should be 8,6,1,7
        intvec edgevec(3) = 7,9,10,8;
        intvec edgevec(4) = 11,10,9,12;
        intvec edgevec(5) = 12,4,5,11;
        intvec edgevec(6) = 5,3,2,6;
        intvec vertexParity = -1,1,1,1,1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }
    
    if( s == "6_2" )
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 4,5,6,1;
        intvec edgevec(3) = 3,7,8,5;
        intvec edgevec(4) = 9,10,7,2;
        intvec edgevec(5) = 11,12,10,9;
        intvec edgevec(6) = 6,8,12,11;
        intvec vertexParity = -1,-1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }
    
    if( s == "6_3" )
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 9,5,1,10;
        intvec edgevec(3) = 7,6,2,5;
        intvec edgevec(4) = 8,3,6,7;
        intvec edgevec(5) = 10,4,12,11;
        intvec edgevec(6) = 11,12,8,9;
        intvec vertexParity = 1,-1,1,1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }
    
    // 7 crossing knots ////////////////////////////////////////
    
    if( s == "7_1" )
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 14,13,2,1;
        intvec edgevec(3) = 12,11,13,14;
        intvec edgevec(4) = 10,9,11,12;
        intvec edgevec(5) = 8,7,9,10;
        intvec edgevec(6) = 6,5,7,8;
        intvec edgevec(7) = 4,3,5,6;
        intvec vertexParity = 1,1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6),edgevec(7));
    }
    
    // 8 crossing knots ////////////////////////////////////////
    
    if( s == "8_19" )
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 5,6,7,2;
        intvec edgevec(3) = 4,3,8,9;
        intvec edgevec(4) = 8,7,10,11;
        intvec edgevec(5) = 12,13,5,1;
        intvec edgevec(6) = 14,15,13,12;
        intvec edgevec(7) = 9,11,16,14;
        intvec edgevec(8) = 16,10,6,15;
        intvec vertexParity = -1,-1,-1,-1,-1,-1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6),edgevec(7),edgevec(8));
    }

    // 11 crossing knots ////////////////////////////////////////
    
    if( s == "11KTC_1" )
    {
        intvec edgevec(1) = 1,4,3,2;
        intvec edgevec(2) = 2,7,6,5;
        intvec edgevec(3) = 3,8,22,7;
        intvec edgevec(4) = 5,6,9,12;
        intvec edgevec(5) = 9,22,10,13;
        intvec edgevec(6) = 12,18,15,14;
        intvec edgevec(7) = 14,15,17,16;
        intvec edgevec(8) = 16,17,21,1;
        intvec edgevec(9) = 13,20,19,18;
        intvec edgevec(10) = 21,19,20,11;
        intvec edgevec(11) = 11,10,8,4;
        intvec vertexParity = 1,1,-1,1,-1,-1,-1,-1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6),edgevec(7),edgevec(8),edgevec(9),edgevec(10),edgevec(11));
    }
    
    // LINKS ////////////////////////////////////////
    
    if( s == "hopf" || s == "L2a1" )
    {
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 4,3,2,1;
        intvec vertexParity = -1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2));
    }

    if( s == "hopf_version_2" || s == "L2a1_version_2" )
    {
        intvec edgevec(1) = 3,2,1,4;
        intvec edgevec(2) = 4,1,2,3;
        intvec vertexParity = -1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2));
    }
    
    if( s == "solomon" || s == "L4a1")
    {    
        intvec edgevec(1) = 1,2,3,4;
        intvec edgevec(2) = 8,3,2,7;
        intvec edgevec(3) = 4,5,6,1;
        intvec edgevec(4) = 7,6,5,8;
        intvec vertexParity = 1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4));
    }

    if( s == "solomon_version_2" || s == "L4a1_version_2")
    {
        intvec edgevec(1) = 2,3,4,1;
        intvec edgevec(2) = 7,8,3,2;
        intvec edgevec(3) = 1,4,5,6;
        intvec edgevec(4) = 6,5,8,7;
        intvec vertexParity = 1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4));
    }

    if( s == "whitehead" || s == "L5a1" )
    {
        intvec edgevec(1) = 4,3,1,2;
        intvec edgevec(2) = 2,1,10,9;
        intvec edgevec(3) = 8,10,3,5;
        intvec edgevec(4) = 6,5,4,7;
        intvec edgevec(5) = 7,9,8,6;
        intvec vertexParity = 1,1,1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5));
    }

    if( s == "whitehead_version_2" || s == "L5a1_version_2" )
    {
        intvec edgevec(1) = 2,4,3,1;
        intvec edgevec(2) = 8,7,5,4;
        intvec edgevec(3) = 9,6,7,8;
        intvec edgevec(4) = 1,10,9,2;
        intvec edgevec(5) = 3,5,6,10;
        intvec vertexParity = -1,1,1,-1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5));
    }

    if( s == "borromeanRings" || s == "L6a4" )
    {    
        intvec edgevec(1) = 1,2,3,7;
        intvec edgevec(2) = 3,4,5,6;
        intvec edgevec(3) = 7,6,8,9;
        intvec edgevec(4) = 11,10,4,2;
        intvec edgevec(5) = 8,5,10,12;
        intvec edgevec(6) = 9,12,11,1;
        intvec vertexParity = -1,1,-1,1,1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "borromeanRings_version_2" || s == "L6a4_version_2" )
    {    
        intvec edgevec(1) = 1,2,3,7;
        intvec edgevec(2) = 6,3,4,5;
        intvec edgevec(3) = 9,7,6,8;
        intvec edgevec(4) = 10,4,2,11;
        intvec edgevec(5) = 8,5,10,12;
        intvec edgevec(6) = 12,11,1,9;
        intvec vertexParity = -1,-1,1,1,-1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "borromeanRings_version_3" || s == "L6a4_version_3" )
    {    
        intvec edgevec(1) = 2,3,7,1;
        intvec edgevec(2) = 5,6,3,4;
        intvec edgevec(3) = 9,7,6,8;
        intvec edgevec(4) = 10,4,2,11;
        intvec edgevec(5) = 12,8,5,10;
        intvec edgevec(6) = 11,1,9,12;
        intvec vertexParity = 1,1,1,-1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "borromeanRings_version_4" || s == "L6a4_version_4" )
    {    
        intvec edgevec(1) = 2,3,7,1;
        intvec edgevec(2) = 4,5,6,3;
        intvec edgevec(3) = 7,6,8,9;
        intvec edgevec(4) = 11,10,4,2;
        intvec edgevec(5) = 12,8,5,10;
        intvec edgevec(6) = 1,9,12,11;
        intvec vertexParity = 1,-1,-1,1,-1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a5")
    {    
        intvec edgevec(1) = 2,1,12,3;
        intvec edgevec(2) = 3,5,4,2;
        intvec edgevec(3) = 11,6,5,10;
        intvec edgevec(4) = 7,4,6,8;
        intvec edgevec(5) = 8,9,1,7;
        intvec edgevec(6) = 10,12,9,11;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a5_version_1") // use x(1) vor mark 1, x(8) vor mark 2, x(7) vor mark 3
    {    
        intvec edgevec(1) = 1,6,12,2;
        intvec edgevec(2) = 2,3,4,1;
        intvec edgevec(3) = 7,4,5,8;
        intvec edgevec(4) = 8,9,6,7;
        intvec edgevec(5) = 11,5,3,10;
        intvec edgevec(6) = 10,12,9,11;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a5_version_2") // use x(1) vor mark 1, x(8) vor mark 2, x(7) vor mark 3
    {    
        intvec edgevec(1) = 2,1,6,12;
        intvec edgevec(2) = 3,4,1,2;
        intvec edgevec(3) = 8,7,4,5;
        intvec edgevec(4) = 9,6,7,8;
        intvec edgevec(5) = 11,5,3,10;
        intvec edgevec(6) = 10,12,9,11;
        intvec vertexParity = -1,-1,-1,-1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a5_version_3") // use x(1) vor mark 1, x(8) vor mark 2, x(7) vor mark 3
    {    
        intvec edgevec(1) = 6,12,2,1;
        intvec edgevec(2) = 1,2,3,4;
        intvec edgevec(3) = 7,4,5,8;
        intvec edgevec(4) = 8,9,6,7;
        intvec edgevec(5) = 5,3,10,11;
        intvec edgevec(6) = 11,10,12,9;
        intvec vertexParity = -1,-1,1,1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a5_version_4") // use x(1) vor mark 1, x(8) vor mark 2, x(7) vor mark 3
    {    
        intvec edgevec(1) = 1,6,12,2;
        intvec edgevec(2) = 2,3,4,1;
        intvec edgevec(3) = 4,5,8,7;
        intvec edgevec(4) = 7,8,9,6;
        intvec edgevec(5) = 10,11,5,3;
        intvec edgevec(6) = 12,9,11,10;
        intvec vertexParity = 1,1,-1,-1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6n1" )
    {    
        intvec edgevec(1) = 4,9,1,3;
        intvec edgevec(2) = 7,8,4,5;
        intvec edgevec(3) = 10,12,9,8;
        intvec edgevec(4) = 5,3,2,6;
        intvec edgevec(5) = 6,11,10,7;
        intvec edgevec(6) = 2,1,12,11;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6n1_version_1" )
    {    
        intvec edgevec(1) = 1,5,12,11;
        intvec edgevec(2) = 2,11,10,9;
        intvec edgevec(3) = 3,4,1,2;
        intvec edgevec(4) = 9,8,6,3;
        intvec edgevec(5) = 6,7,5,4;
        intvec edgevec(6) = 10,12,7,8;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6n1_version_2" )
    {    
        intvec edgevec(1) = 5,12,11,1;
        intvec edgevec(2) = 2,11,10,9;
        intvec edgevec(3) = 4,1,2,3;
        intvec edgevec(4) = 3,9,8,6;
        intvec edgevec(5) = 6,7,5,4;
        intvec edgevec(6) = 8,10,12,7;
        intvec vertexParity = -1,1,-1,-1,1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6n1_version_3" )
    {    
        intvec edgevec(1) = 1,5,12,11;
        intvec edgevec(2) = 11,10,9,2;
        intvec edgevec(3) = 2,3,4,1;
        intvec edgevec(4) = 9,8,6,3;
        intvec edgevec(5) = 4,6,7,5;
        intvec edgevec(6) = 12,7,8,10;
        intvec vertexParity = 1,-1,-1,1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6n1_version_4" )
    {    
        intvec edgevec(1) = 11,1,5,12;
        intvec edgevec(2) = 9,2,11,10;
        intvec edgevec(3) = 3,4,1,2;
        intvec edgevec(4) = 8,6,3,9;
        intvec edgevec(5) = 7,5,4,6;
        intvec edgevec(6) = 10,12,7,8;
        intvec vertexParity = -1,-1,1,-1,-1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a1" )
    {    
        intvec edgevec(1) = 3,1,5,4;
        intvec edgevec(2) = 4,6,2,3;
        intvec edgevec(3) = 2,11,12,1;
        intvec edgevec(4) = 10,12,11,9;
        intvec edgevec(5) = 8,9,6,7;
        intvec edgevec(6) = 7,5,10,8;
        intvec vertexParity = 1,1,-1,-1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a1_version_2" )
    {    
        intvec edgevec(1) = 1,5,4,3;
        intvec edgevec(2) = 3,4,6,2;
        intvec edgevec(3) = 2,11,12,1;
        intvec edgevec(4) = 10,12,11,9;
        intvec edgevec(5) = 9,6,7,8;
        intvec edgevec(6) = 8,7,5,10;
        intvec vertexParity = -1,-1,-1,-1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a2" )
    {    
        intvec edgevec(1) = 4,2,1,3;
        intvec edgevec(2) = 11,1,2,5;
        intvec edgevec(3) = 3,12,6,4;
        intvec edgevec(4) = 5,6,8,7;
        intvec edgevec(5) = 7,8,10,9;
        intvec edgevec(6) = 9,10,12,11;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a2_version_2" )
    {    
        intvec edgevec(1) = 3,4,2,1;
        intvec edgevec(2) = 1,2,5,11;
        intvec edgevec(3) = 12,6,4,3;
        intvec edgevec(4) = 7,5,6,8;
        intvec edgevec(5) = 8,10,9,7;
        intvec edgevec(6) = 11,9,10,12;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a3" )
    {    
        intvec edgevec(1) = 1,2,12,11;
        intvec edgevec(2) = 3,4,2,1;
        intvec edgevec(3) = 5,6,4,3;
        intvec edgevec(4) = 7,8,6,5;
        intvec edgevec(5) = 9,10,8,7;
        intvec edgevec(6) = 11,12,10,9;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L6a3_version_2" )
    {    
        intvec edgevec(1) = 11,1,2,12;
        intvec edgevec(2) = 4,2,1,3;
        intvec edgevec(3) = 3,5,6,4;
        intvec edgevec(4) = 8,6,5,7;
        intvec edgevec(5) = 7,9,10,8;
        intvec edgevec(6) = 12,10,9,11;
        intvec vertexParity = 1,1,1,1,1,1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6));
    }

    if( s == "L7a7" )
    {    
        intvec edgevec(1) = 5,1,7,6;
        intvec edgevec(2) = 6,8,4,5;
        intvec edgevec(3) = 3,11,1,2;
        intvec edgevec(4) = 2,4,12,3;
        intvec edgevec(5) = 12,13,14,11;
        intvec edgevec(6) = 9,10,13,8;
        intvec edgevec(7) = 7,14,10,9;
        intvec vertexParity = 1,1,1,1,-1,-1,-1;
        retLink = vertexParity, list(edgevec(1),edgevec(2),edgevec(3),edgevec(4),edgevec(5),edgevec(6),edgevec(7));
    }

    return(retLink);
}

////////////////////////////////////////////////////////////////////
// linkNumberOfVariables
//
// Given a link (as returned by linkLibrary) returns the number of
// variables needed in the basering

proc linkNumberOfVariables(list linkData)
{
    list edgeList = linkData[2];
    
    int i;
    int mx = 1;
    for(i=1;i<=size(edgeList);i++)
    {
        int j;
        for(j=1;j<=4;j++)
        {
            if( edgeList[i][j] > mx )
            {
                mx = edgeList[i][j];
            }
        }
    }
    
    return(mx);
}

////////////////////////////////////////////////////////////////////
// linkParity
//
// Given a link (as returned by linkLibrary) returns the parity, in
// the sense of [KR]. Namely, we replace all crossings with the smooth
// resolution and count the number of circles in the resulting link.
// The "parity" is this integer modulo 2. i.e. we return 0 or 1.

proc linkParity(list linkData)
{
    // We maintain a list of "unclaimed" edges, and cross off
    // from this list as we enumerate components using
    // enumerateComponentFromEdge
    
    list unclaimedEdges;
    int numVar = linkNumberOfVariables(linkData);
    
    int i;
    for(i=1; i<=numVar; i++)
    {
        unclaimedEdges = unclaimedEdges + list(i);
    }
    
    list componentList;
    
    while( size(unclaimedEdges) != 0 )
    {
        list thisComp = enumerateComponentFromEdgeInSmoothing( linkData, unclaimedEdges[1] );
        
        componentList = componentList + list(thisComp);
        
        // Now delete from our unclaimed list all the edges
        // we found in this component
        int i;
        for(i=1;i<=size(thisComp);i++)
        {
            int j;
            for(j=1;j<=size(unclaimedEdges);j++)
            {
                if( unclaimedEdges[j] == thisComp[i] )
                {
                    unclaimedEdges = delete(unclaimedEdges, j);
                    break;
                }
            }
        }
        
        kill i,j,thisComp;
    }
    
    int parity = size(componentList) mod 2;
    return(parity);
}

////////////////////////////////////////////////////////////////////
// enumerateComponentFromEdge
//
// Given a link (as returned by linkLibrary) and an edge in the
// link, returns the list of edges in the same component of the
// link as the given edge. The ordering follows the orientation.

proc enumerateComponentFromEdge(list linkData, int m)
{
    list comp;
    int currentEdge = m;
    
    while( 1 )
    {
        comp = comp + list(currentEdge);
        currentEdge = nextEdgeFromEdge(linkData, currentEdge);
        
        if( currentEdge == m ){ break; }
    }
    
    return(comp);
}

////////////////////////////////////////////////////////////////////
// enumerateComponentFromEdgeInSmoothing
//
// Given a link (as returned by linkLibrary) and an edge in the
// link, returns the list of edges in the same component of the
// "smoothed" link obtained by replacing every crossing by its
// 0-smoothing in the sense of resolutions.

proc enumerateComponentFromEdgeInSmoothing(list linkData, int m)
{
    list comp;
    int currentEdge = m;
    
    while( 1 )
    {
        comp = comp + list(currentEdge);
        currentEdge = nextEdgeFromEdgeInSmoothing(linkData, currentEdge);
        
        if( currentEdge == m ){ break; }
    }
    
    return(comp);
}

////////////////////////////////////////////////////////////////////
// nextEdgeFromEdge
//
// Given a link (as returned by linkLibrary) and an edge in the
// link, returns the next edge following the orientation of the link.

proc nextEdgeFromEdge(list linkData, int m)
{
    list edges = linkData[2];
    
    int i;
    for(i=1; i<=size(edges); i++)
    {
        intvec e = edges[i];
        
        if( e[3] == m )
        {
            // Our edge is the bottom right leg, so return the top left leg
            return(e[1]);
        }
        
        if( e[4] == m )
        {
            // Our edge is the bottom left leg, so return the top right leg
            return(e[2]);
        }
        
        kill e;
    }
    
    print("[nextEdgeFromEdge] Error: could not find next edge.");
    return(-1);
}

////////////////////////////////////////////////////////////////////
// nextEdgeFromEdgeInSmoothing
//
// Given a link (as returned by linkLibrary) and an edge in the
// link, returns the next edge if we erase the crossing and replace
// it with the smoothing (i.e. the 0 state).

proc nextEdgeFromEdgeInSmoothing(list linkData, int m)
{
    list edges = linkData[2];
    
    int i;
    for(i=1; i<=size(edges); i++)
    {
        intvec e = edges[i];
        
        if( e[3] == m )
        {
            // Our edge is the bottom right leg, so return the top right leg
            return(e[2]);
        }
        
        if( e[4] == m )
        {
            // Our edge is the bottom left leg, so return the top left leg
            return(e[1]);
        }
        
        kill e;
    }
    
    print("[nextEdgeFromEdgeInSmoothing] Error: could not find next edge.");
    return(-1);
}

/////////////////////////////////////////////////////////////////////////////////////
// gpoly - returns the polynomial g of [KR] p.48.
/////////////////////////////////////////////////////////////////////////////////////

proc gpoly(poly s1, poly s2, int n)
{
    poly g = s1^(n+1);
    
    int i;
    for(i=1;i<=(n+1) div 2;i++)
    {
        g = g + (n+1) * (-1)^(i) * ( binomial(n-i,i-1) * s2^(i) * s1^(n+1-2*i) )/i;
    }
    
    return(g);
}

////////////////////////////////////////////////////////////////////
// mfWideEdge
//
// We are given four variables a,b,c,d and an integer N, and we return
// the matrix factorisation of a^(n+1) + b^(n+1) - c^(n+1) - d^(n+1)
// associated by Khovanov and Rozansky to a wide edge in the planar
// diagram of a link (see p.48 of [KR]). Let us include the
// following picture of a wide edge
//
//       a        b
//         \   /
//          \/
//          |
//          |
//         /\
//       /   \
//     d       c
//
// where a,b,c,d are variables labelled on the edges as shown. Note that this
// is the ordering convention of Figure 8 on p.9 of [KR]. With this convention
// if you call
//
// matrix D = mfWideEdge(x(1),x(2),x(3),x(4),N);
//
// Then the matrix D you get will be the (Q1,Q2) pair written on p.50 of [KR]. Note however
// that everything is symmetric for the exchange x(1) <-> x(2) and x(4) <-> x(3).
////////////////////////////////////////////////////////////////////

proc mfWideEdge(poly a, poly b, poly c, poly d, int n)
{
    poly W = a^(n+1) + b^(n+1) - c^(n+1) - d^(n+1);
    
    // Following KR we write W = (a + b - c - d)u1 + (ab - cd)u2:
    poly u1 = (gpoly(a + b, a*b, n) - gpoly(c + d, a*b, n))/(a+b-c-d);
    poly u2 = (gpoly(c + d, a*b, n) - gpoly(c + d, c*d, n))/(a*b-c*d);
    
    matrix koszul1[2][2] = 0, a + b - c - d, u1, 0;
    matrix koszul2[2][2] = 0, a * b - c * d, u2, 0;
    
    // Note the order that we tensor these two Koszul complexes is opposite
    // to the order suggested on p.48 of [KR]. This is done so as to agree
    // with the matrices written on p.50.
    matrix D = MFtensor(koszul2, koszul1);
    
    return(D);
}

proc mfWideEdgeGrading(int n)
{
    // Note that deg(u1) = n and deg(u2) = n - 1 so as written on 
    // [KR, p.48] the Koszul factorisations koszul1 and koszul2 have
    // respective grading vectors
    intvec g = 0, 1-n;
    intvec h = 0, 3-n;
    
    // Grading vector on the tensor product:
    intvec r = MFtensorGradings(h, g);
    
    // The grading is then shifted by -1, as defined on [KR, p.48]:
    int i;
    for(i=1;i<=size(r);i++){ r[i] = r[i] - 1; }
    
    return(r);
}

////////////////////////////////////////////////////////////////////
// mfIdentityDefect
//
// We are given two potentials W1 and W2 and two variables y and x.
// We assume that W1 = W(y) and W2 = W(x) for some fixed polynomial
// W and we return the identity defect W2 -> W1.
////////////////////////////////////////////////////////////////////

proc mfIdentityDefect(poly W1, poly W2, poly y, poly x)
{
    matrix koszul1[2][2] = 0, y - x, (W1 - W2)/(y - x), 0;
    return(koszul1);
}

// We are given the degree of some homogeneous potential.
proc mfIdentityDefectGrading(int dg)
{
    int n = dg - 1;
    
    // We have |W| = 2(n+1) and the identity defect is given the grading
    // R ---> R{1+n-|(W1-W2)/(y-x)|} --> R, but |(W1-W2)/(y-x)| = 2 * n so
    // this is R ---> R{1-n} ---> R, see [KR, p.45].
    
    intvec g = 0, 1-n;
    return(g);
}

////////////////////////////////////////////////////////////////////
// mfIdentityEdge
//
// We are given four variables a,b,c,d and an integer N, and we return
// the matrix factorisation of a^(n+1) + b^(n+1) - c^(n+1) - d^(n+1)
// associated by Khovanov and Rozansky to the Gamma0 diagram of Figure 8.
// Let us include the following picture of a wide edge
//
//       a        b
//         \   /
//          \/
//          |
//          |
//         /\
//       /   \
//     d       c
//
// where a,b,c,d are variables labelled on the edges as shown. The return
// is the tensor product of L^b_c with L^a_d (in KR's notation), in that
// order. Note that this is not symmetric with respect to the interchange
// a <-> b or c <-> d, so the ordering of the variables matters.
////////////////////////////////////////////////////////////////////

proc mfIdentityEdge(poly a, poly b, poly c, poly d, int n)
{
    matrix Aa = mfIdentityDefect( a^(n+1), d^(n+1), a, d );
    matrix Ab = mfIdentityDefect( b^(n+1), c^(n+1), b, c );
    
    // Note the order of the tensor factors (this is opposite to what
    // a first reading of the third line of [KR, p.50] would suggest,
    // but is necessary to reproduce the matrices P0 and P1 there, and
    // to make chi0 = (U0,U1) into a morphism of MFs).
    matrix A = MFtensor(Ab, Aa);
    
    return(A);
}

proc mfIdentityEdgeGrading(int n)
{
    intvec Aagr = mfIdentityDefectGrading( n+1 );
    intvec Abgr = mfIdentityDefectGrading( n+1 );
    intvec Agr = MFtensorGradings(Abgr, Aagr);
    
    return(Agr);
}

////////////////////////////////////////////////////////////////////
// checklist
//
// checklist(l,p) gives 1 if p is an element of the list l,
// otherwise it gives 0.
////////////////////////////////////////////////////////////////////

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


////////////////////////////////////////////////////////////////////
// removefromlist
//
// removefromlist(l,p) gives back the list l, but with all instances of
// entries p deleted. For example:
// removefromlist( list(1,2,2,3,2), 2 ) = list(1,3). 
////////////////////////////////////////////////////////////////////

proc removefromlist(list l, poly p)
{
    for(int i=1; i<=size(l); i++)
    {
        def li = l[i];
        if(p==li)
        {
            l = delete(l,i);
            i=i-1;
        }
    }
    return(l);
}


////////////////////////////////////////////////////////////////////
// mfVertex
//
// mfVertex(zeroOrOne, varis, N) gives the "KR-vertex" of type 0
// (tensor product of two identity defects) or of type 1 (wide edge
// defect) in the list of variables varis. N as in sl(N). The ordering
// of the variables in varis corresponds to 
//
//     varis[1]   varis[2]
//           \   /
//            \/
//            |
//            |
//           /\
//         /   \
//   varis[4]    varis[3]
//
// See the comments preceeding mfIdentityEdge for an idea of why the
// ordering here is significant.
////////////////////////////////////////////////////////////////////

proc mfVertex(int zeroOrOne, list varis, int N)
{
    if( size(varis) != 4 )
    {
        print("[mfVertex] Passed a variable list of incorrect size");
        return();
    }
    
    if( zeroOrOne == 0 )
    {
        matrix A = mfIdentityEdge( varis[1], varis[2], varis[3], varis[4], N );
        return(A);
    }
        
    if( zeroOrOne == 1 )
    {
        matrix A = mfWideEdge( varis[1], varis[2], varis[3], varis[4], N );
        return(A);
    }
}

////////////////////////////////////////////////////////////////////
// mfVertexGrading
////////////////////////////////////////////////////////////////////

proc mfVertexGrading(int zeroOrOne, int N)
{
   if( zeroOrOne == 0 )
   {
       intvec Agr = mfIdentityEdgeGrading(N);
       return(Agr);
   }
   
   if( zeroOrOne == 1 )
   {
       intvec Agr = mfWideEdgeGrading(N);
       return(Agr);
   }
}

////////////////////////////////////////////////////////////////////
// mfChi
//
// mfChi(zeroOrOne,p(1),p(2),p(3),p(4),N) is the map C(Gamma^0) --> C(Gamma^1)
// from [KR, p.50] for zeroOrOne=0, and the map C(Gamma^1) --> C(Gamma^0)
// from [KR, p.51] for zeroOrOne=1.
//
// NOTE: The definition of the maps in [KR] involves params mu, lambda,
// which are allowed to take any integer value. We set these equal to zero,
// but any other value gives a morphism too. 
////////////////////////////////////////////////////////////////////

proc mfChi(int zeroOrOne, poly p(1), p(2), p(3), p(4), int N)
{
    if( zeroOrOne == 0 )
    {
        number mu = 0;

        poly pi23 = ( p(2)^(N+1) - p(3)^(N+1) )/( p(2)-p(3) );
        poly u1 = ( gpoly(p(1)+p(2), p(1)*p(2), N) - gpoly(p(3)+p(4), p(1)*p(2), N) )/( p(1)+p(2)-p(3)-p(4) );
        poly u2 = ( gpoly(p(3)+p(4), p(1)*p(2), N) - gpoly(p(3)+p(4), p(3)*p(4), N) )/( p(1)*p(2)-p(3)*p(4) );
        poly aa1 = (mu-1) * u2 + ( u1 +  p(1) * u2 - pi23 )/( p(1) - p(4) );        
        
        matrix chi0[4][4] = p(4)-p(2) + mu*(p(1)+p(2)-p(3)-p(4)), 0, 0,    0,
                        aa1,        1, 0,    0,
                        0,         0, p(4) + mu*(p(1)-p(4)), -p(2) + mu*(p(2)-p(3)),
                        0,         0, -1,   1;
        return(chi0);
    }

    if( zeroOrOne == 1 )
    {
        number lambda = 0;
        
        poly pi23 = ( p(2)^(N+1) - p(3)^(N+1) )/( p(2)-p(3) );
        poly u1 = ( gpoly(p(1)+p(2), p(1)*p(2), N) - gpoly(p(3)+p(4), p(1)*p(2), N) )/( p(1)+p(2)-p(3)-p(4) );
        poly u2 = ( gpoly(p(3)+p(4), p(1)*p(2), N) - gpoly(p(3)+p(4), p(3)*p(4), N) )/( p(1)*p(2)-p(3)*p(4) );
        poly aa2 = lambda*u2 + ( u1 + p(1) * u2 - pi23 )/( p(4)-p(1) );
        poly aa3 = lambda*(p(3)+p(4)-p(1)-p(2)) + p(1) - p(3);

        matrix chi1[4][4] = 1, 0, 0,    0,
                        aa2,        aa3, 0,    0,
                        0,         0, 1, p(3) + lambda*(p(2)-p(3)),
                        0,         0, 1, p(1) + lambda*(p(4)-p(1));
    
        return(chi1);
    }
}

/////////////////////////////////////////
// powervec
//
// Given an integer m returns a list of all intvecs of length m with entries 0,1.

proc powervec(int m)
{
    if( m == 1 )
    {
        intvec a = 0;
        intvec b = 1;
        list r = a,b;
        kill a,b;
        return(r);
    }
    
    list l = powervec(m-1);
    list r;
    
    int i;
    for(i=1;i<=size(l);i++)
    {
       // Prepend either 0 or 1 to each entry of l:
        intvec v = l[i];
        intvec a = 0,v;
        intvec b = 1,v;
        r = r + list(a,b);
        kill v,a,b;
    }
    
    kill l;
    return(r):
}

/////////////////////////////////////////
// intvecToString
//
// NOTE that in mfCompileLink below we only need to evaluate intvecToString on
// intvecs that are made up only of 0s and 1s, so there cannot be a problem with
// the fact that intvecToString maps e.g. both 11,2,3 and 1,2,3,4 to the same string.
/////////////////////////////////////////

proc intvecToString(intvec v)
{
    string s;
    int i;
    for( i = 1; i <= size(v); i++ )
    {
        s = s + string(v[i]);
    }
    
    return(s);
}

proc externalDegreeOfWeb(intvec v, intvec vertexParity)
{
    // 0 is always degree zero, and 1 is either degree -1 or +1 depending on the crossing type.
    int i;
    int d;
    for( i = 1; i <= size(v); i++ )
    {
        d = d + v[i] * vertexParity[i];
    }
    
    return(d);
}

////////////////////////////////////////////////////////////////////
// superblockmat
//
// We are given a list of matrices which are supposed to fit into
// a block matrix with "numCols" columns, and we return the block matrix.
// The list is interpreted as a block matrix as follows: the first numCols
// entries form the first row, the second numCols entries the second row, etc.

proc superblockmat(int numCols, list #)
{
    if( size(#) mod numCols != 0 )
    {
        print("[superblockmat] The number of given matrices is not divisble by the number of columns, exiting.");
        return();
    }
    
    int numRows = size(#) div numCols;
    
    // We first collect each row of the block matrix into a list.
    list l;
    
    int i,j,k,t,s;
    for(s = 1;s<=numRows;s++)
    {
        // In this row of the block matrix we have the matrices
        // #[(s-1)*numCols + 1], ...., #[s * numCols].
        int b = nrows(#[s*numCols]);
        
        for(j=1;j<=b;j++)
        {
            list thisrow;
            
            for(k=1;k<=numCols;k++)
            {
                // We are now inside #[(s-1)*numCols + k].
                matrix cMat = #[(s-1)*numCols +k];
                for(t=1;t<=ncols(cMat);t++)
                {
                    thisrow = thisrow + list(cMat[j,t]);
                }
                kill cMat;
            }
            
            l = l + list(thisrow);
        }
    }
    
    // Convert l to a matrix:
    matrix A[size(l)][size(l[1])];
    for(i=1;i<=size(l);i++)
    {
        for(j=1;j<=size(l[1]);j++)
        {
            A[i,j] = l[i][j];
        }
    }
    
    return(A);
}

////////////////////////////////////////////////////////////////////
// computeCohomology
//
// Takes as input a complex of graded vector spaces, specified as a
// tuple (m,L,G) as described as the ouput of KRcomplex. The return is
// a tuple (m,V,H) where V is a list giving the dimension of the
// cohomologies in each degree, and H is a list giving the weights of
// each of these cohomology spaces.

proc computeCohomology(list cpx)
{
    list C = cpx[2];
    int minDegree = cpx[1];
    int maxDegree = minDegree + size(C);
    list weightsPerCohomologicalDegree = cp[3];
    
    def cohom(1) = rightKernel(C[1]);
    int i;
    for( i=2; i<=maxDegree-minDegree; i++ )
    {
        def ker(i) = rightKernel(C[i]);
        def im(i-1) = module(C[i-1]);
        module coho(i) = reduce(ker(i),std(im(i-1)));
        def cohom(i) = std( coho(i) );
    }
    def ker(maxDegree-minDegree+1) = module( unitmat(size(weightsPerCohomologicalDegree[maxDegree-minDegree+1])));
    def im(maxDegree-minDegree) = module(C[maxDegree-minDegree]);
    module coho(maxDegree-minDegree+1) = reduce(ker(maxDegree-minDegree+1),std(im(maxDegree-minDegree)));
    def cohom(maxDegree-minDegree+1) = std( coho(maxDegree-minDegree+1) );

    list V;
    list H;
    int j, k;
    for( i=1; i<=maxDegree-minDegree+1; i++)
    {
        V[i] = size(cohom(i));
        list wt;
        
        for( j=1; j<=size(cohom(i)); j++ )
        {
            vector v = cohom(i)[j];
            
            if( v == [0] )
            {
                print("[computeCohomology] Zero generator in cohomology, exiting.");
                return();
            }
            
            int vwt;
            int found = 0;
            
            // Look at which generators occur with nonzero coefficients
            // in v, and from this deduce deg(v). Note that each generator
            // occurring with a nonzero coeff should have the same degree,
            // otherwise v has no well-defined degree.
            for( k=1; k<=nrows(v); k++ )
            {
                if( v[k] != 0 )
                {
                    if( found == 0 )
                    {
                        found = 1;
                        vwt = weightsPerCohomologicalDegree[i][k];
                    }
                    else
                    {
                        if( vwt != weightsPerCohomologicalDegree[i][k] )
                        {
                            print("[computeCohomology] Generator " + string(j) +
                                " of cohom(" + string(i) + ") has no well-defined degree, exiting.");
                            print("Generator: " + string(v) );
                            print("Grading:   " + string(weightsPerCohomologicalDegree[i]));
                            return();
                        }
                    }
                }
            }
            
            // If we are still here, we found a consistent degree.
            wt[j] = vwt;
        }
        
        H[i] = wt;
    }

    return(list(minDegree,V,H));
}

////////////////////////////////////////////////////////////////////
// gradedPoincare
//
// Takes as input a sequence of graded vector spaces, specified as a
// tuple (m,V,H) as described as the ouput of computeCohomology. The return
// is the graded Poincare polynomial. The second and third argument are 
// the params to use (i.e., the q and t).
//
// Recall that given a bigraded vector space H, the Poincare polynomial is
//
//     sum_i,j t^i * q^j * dim( H^(i,j) )
//
// where i, j range over all integers, i is a cohomological degree and j
// is the weight or internal degree.

proc gradedPoincare(list cohom, number v, number u)
{
    int minDegree = cohom[1];
    list V = cohom[2]; // Dimensions of cohomologies
    list H = cohom[3]; // Weight vectors of cohomologies

    number p;

    int i,j;
    for( i=1; i<=size(V); i++)
    {
        for( j=1; j<=size(H[i]); j++ )
        {
            p = p + u^(i + minDegree - 1) * v^(H[i][j]);
        }
    }

    return(p);
}

////////////////////////////////////////////////////////////////////
// gradedEuler
//
// Takes as input a sequence of graded vector spaces, specified as a
// tuple (m,V,H) as described as the ouput of computeCohomology. The return
// is the graded Euler characteristic. The second argument is the param
// to use (i.e., the q).

proc gradedEuler(list cohom, number v)
{
    int minDegree = cohom[1];
    list V = cohom[2]; // Dimensions of cohomologies
    list H = cohom[3]; // Weight vectors of cohomologies
    
    // Now we want to compute the Euler char as
    //     sum_i,j (-1)^i * q^j * dim( H^(i,j) )
    // where is the cohomological degree, and j is the weight.

    number euler;

    int i,j;
    for( i=1; i<=size(V); i++)
    {
        for( j=1; j<=size(H[i]); j++ )
        {
            int sign = (i + minDegree + 1) mod 2;
            euler = euler + (-1)^(sign) * v^(H[i][j]);
        }
    }

    return(euler);
}

////////////////////////////////////////////////////////////////////
// KRcomplex
// 
// Given a link, defined in the format specified prior to mfCompileLink,
// and the "compilation" h of this link as given by mfCompileLink, we
// return a pair (m,L,G) where m is the minimum nonzero degree of the KR
// complex, L is the list of matrices giving the differentials, and G is
// the list of intvecs giving the weights on each of the vector spaces.

proc KRcomplex(list linkData, list linkComp, int N)
{
    // linkComp[1] is a list of tuples (intvec state, int rank, int degree, intvec grading)
    // limkComp[2] is a list of tuples (intvec state1,intvec state2,int chitype,int chipos,matrix F).
    
    list h = linkComp[2];
    intvec vertexParity = linkData[1];
    
    // Discover which webs in the KR complex live in which cohomological degrees:
    list r = listWebsByDegree(linkData, linkComp, N);
    
    int minDegree = r[1];
    list websByDegree = r[2];
    int maxDegree = minDegree + size(websByDegree) - 1;
    
    // Find the weights (="internal degrees") of the individual, total pieces in each
    // cohomological degree. 
    list weightsPerCohomologicalDegree;
    int i, j;
    for( j=minDegree; j<=maxDegree; j++ )
    {
        int sj = j - minDegree + 1;
        intvec temp = websByDegree[sj][1][3];

        for( i = 2; i<=size(websByDegree[sj]); i++ )
        {
            list currWeb = websByDegree[sj][i];
            temp = temp,currWeb[3];
            kill currWeb;
        }

        weightsPerCohomologicalDegree[sj] = temp;
        kill temp;
    }

    // Try computing the differentials in the KR complex
    list KRKomplex;
    int c;
    for( c=minDegree; c<=maxDegree-1; c++ )
    {
        int sc = c - minDegree + 1; // shifted degree
        int i,j;
        list webmaps;

        for(j=1;j<=size(websByDegree[sc+1]);j++)
        {
            for(i=1;i<=size(websByDegree[sc]);i++)
            {
                list source = websByDegree[sc][i];
                list target = websByDegree[sc+1][j];

                // Both source, target have the format (intvec state, int rank, intvec weight)

                // Look for the map between these two webs in h, if it exists:
                matrix A;
                int s;
                int found = 0;
                for(s=1;s<=size(h);s++)
                {
                    if( h[s][1] == source[1] && h[s][2] == target[1] )
                    {
                        // Calculate the Koszul sign:
                        int signf;
                        int chiposition = h[s][4];
                        int b;
                        for(b=1;b<=chiposition-1;b++)
                        {
                            signf = signf + h[s][1][b];
                        }

                        A = (-1)^(signf) * h[s][5];
                        found = 1;
                        
                        // There can only be one map between source and target. Hence:
                        break;
                    }
                }
                
                if( !found )
                {
                    A = matrix(0,target[2],source[2]);
                }

                webmaps = webmaps + list(A);
            }
        }

        matrix final = superblockmat(size(websByDegree[sc]),webmaps);
        KRKomplex = KRKomplex + list(final);
    }

    // Sanity check:
    for( i=minDegree; i<=maxDegree-2; i++ )
    {
        matrix zerom[nrows(KRKomplex[i+1-minDegree+1])][ncols(KRKomplex[i-minDegree+1])];
        if( KRKomplex[i+1-minDegree+1] * KRKomplex[i-minDegree+1] != zerom )
        {
            print("[KRlinkInvariant] KRcomplex is not a complex at cohomological degree " + string(i) + ", exiting");
            return();
        }
    }
    
    return(list(minDegree,KRKomplex,weightsPerCohomologicalDegree));
}

////////////////////////////////////////////////////////////////////
// listWebsByDegree
//
// Given a link, and the "compilation" of this link as given by linkCompile, we
// return a list L with the following format. Let minDegree be the smallest
// degree with a nonzero entry in the KR complex. Then
//
//      L[i] = list of all triples (A,B,C) where A is a state diagram intvec
//      with "external" or cohomological degree minDegree + i - 1, B
//      is the dimension of the vector space which is the compilation of the
//      corresponding web, and C is an intvec that gives the weights of this vector space.
//
// NOTE: The calculation of the weight is a little subtle, and worth paying attention to.
// The output of mfCompileLink has webs whose weights are just the "bare" ones output by
// the splitting and reduction procedure. But for the KR complex there are additional
// grading shifts defined in [KR,p.81], and we incoporate these shifts here.
//
// NOTE: Webs are ordered in whatever ordering the routine powervec puts them in.
//
// NOTE: We ignore webs that compile to zero dimensional vector spaces.

proc listWebsByDegree(list linkData, list linkComp, int N)
{
    list webComp = linkComp[1];
    list h = linkComp[2];
    
    // webComp is a list of tuples (intvec state, int rank, int degree, intvec grading).
    // h is a list of tuples (intvec state1,intvec state2,int chitype,int chipos,matrix F).
    
    int minDegree = 0;
    int maxDegree = 0;
    
    intvec vertexParity = linkData[1];
    
    int webSiz = size(vertexParity); // number of vertices in our web
    list stateDiagrams = powervec(webSiz);

    int i;
    for( i=1; i<=size(stateDiagrams); i++ )
    {
        intvec web = stateDiagrams[i];
        int d = externalDegreeOfWeb(web, vertexParity);
        if( d < minDegree ){ minDegree = d; }
        if( d > maxDegree ){ maxDegree = d; }
    }

    list websByDegree;
    list e;
    int j;
    
    // TODO: What happens if the complex is all zero in some degree, so that
    // websByDegree[i] is the empty list?
    for( j=minDegree; j<=maxDegree; j++ )
    {
        websByDegree[j-minDegree+1] = e;

        for(i = 1; i <= size(stateDiagrams); i++)
        {
            intvec web = stateDiagrams[i];
            int d = externalDegreeOfWeb(web, vertexParity);
            
            if( d != j ){ i++; continue; }
            
            int rankk = webComp[i][2];

            // We do not include webs compiling to zero.
            if( rankk == 0 ){ print("[listWebsByDegree] Skipping zero web"); i++; continue; }

            // Look for this web as the source or target of something in h:
            intvec weights = webComp[i][4];
            intvec inWeight = weights;
            int c;
            for(c = 1; c <= size(h); c++)
            {
                // Every time we see this web as the source of a map.
                if( h[c][1] == web )
                {    
                    // By [KR,p.81] we have an additional 1 - N or N if the chitype is 0 or 1 resp.
                    int chitype = h[c][3];
            
                    if( chitype == 0 )
                    {
                        weights = weights + 1 - N;
                    }
                    else
                    {
                        weights = weights + N;
                    }
                }
                    
                // Every time we see this web as the target of a map.
                if( h[c][2] == web )
                {
                    // By [KR,p.81] we have an additional -N or N - 1 if the chitype is 0 or 1 resp.
                    int chitype = h[c][3];
            
                    if( chitype == 0 )
                    {
                        weights = weights - N;
                    }
                    else
                    {
                        weights = weights + N - 1;
                    }
                }
            }
            
            websByDegree[j-minDegree+1] = websByDegree[j-minDegree+1] + list(list(web,rankk,weights));
        }
    }
    
    return(list(minDegree,websByDegree));
}

////////////////////////////////////////////////////////////////////
// linkCompile
//
// Calculates the bigraded unreduced sl(N) Khovanov-Rozansky link homology.
//
// Takes as arguments a link, defined as a pair in the sense explained above, and the name
// of a file to use for storing intermediate results. The return is a pair (X,Y) of lists.
// Suppose we have m crossings in our link. A "state diagram" is an assignment of integers 0,1
// to crossings, represented by an intvec with entries 0,1. With this terminology:
//
//   1. X is a list of tuples (intvec state, int rank, int degree, intvec grading) where
//   state ranges over all state diagrams, rank is the dimension of the compilation of the
//   corresponding web, degree is the cohomological degree in which the compilation sits,
//   and grading is the grading vector. NOTE: if we are zero in both degrees then we set
//   degree to -1.
//
//   2. Y is a list of tuples (intvec state1,intvec state2,int chitype,int chipos,matrix F).
//   where state1, state2 range over all state diagrams, chitype is 0 or 1 resp. if the map
//   is a chi0 or chi1 map, and chipos gives the position in the state diagram where the
//   map is inserted. Finally, F is the matrix of the map, the result of calling
//   webCompileMorphism.
//
// To facilitate the computations we store the results of mfCompileWeb to disk.
// To illustrate the scheme, take the state 010 and state 110, let us say that the underlying
// link has a positive crossing in the first position so that there is a map 010 -> 110 of
// chitype 0. Internally to mfCompileLink we would then run
//
//      list map010to110 = mfCompileWeb( suitable params );
//
// This list is then stored to the given file using MPfile.
//
// Before computing a given tuple we load in the contents of the datafile and see if it has
// already been computed.
//
// NOTE: The maps returned by webCompileMorphism are degree one, and we do not mess with this;
// we leave all the grading vectors unchanged.

proc linkCompile(list linkData, string fileName, int N)
{
    // We do not change rings: we assume that the ambient ring has the same number of variables
    // as edges in the link, and we associate ith variable var(i) with the edge i. Of course
    // if there is a save file, then loading this save file will set the ring.
    int useSanityChecks = blowFlags("sanity_checks"); 
    int useMPfile = blowFlags("use_mpfile");
    
    system("--ticks-per-sec",1000);
    int timeLink = timer;
 
    def RRR = basering;
    
    intvec vertexParity = linkData[1];
    list edgeList = linkData[2];
    int webSiz = size(vertexParity);
    
    // Read in the save file, if it exists
    if( status(fileName,"exists") == "yes" )
    {
        if( useMPfile )
        {
            def nR = loadSaveFile(fileName);
            setring nR;
        }
        else
        {
            execute(loadSaveFile(fileName));
        }
    }
        
    list returnListWebs;
    list returnListMaps;
    
    // The compilation strategy is the same for all webs. We always enumerate the variables
    // in the default order, but we have to choose an ordering of the vertices which gives
    // a valid compilation strategy. That is, after fusing the first i vertices, the result
    // should be connected to the (i+1)st vertex by an edge.
    
    // First build the default compilation strategy
    list compStratVert;
    list compStratVar;
    int j;
    for(j=1;j<=size(edgeList); j++)
    {
        compStratVert[j] = j;
    }
    for(j=1;j<=nvars(basering);j++)
    {
        compStratVar = compStratVar + list(var(j));
    }
    list compStrat = compStratVert,compStratVar;
    
    // Now make sure it is valid, and permute the ordering of the vertices if not
    // Note it doesn't matter what MFs we put in the web: compilationStratVerify only cares
    // about the edges
    list stateDiagrams = powervec(webSiz);
    list web = webFromLink(linkData, stateDiagrams[1], N);
    
    if( !compilationStratVerify(web, compStrat) )
    {
        dbprint(printlevel, "[linkCompile] Default compilation strategy is invalid, trying to permute vertices...");
        
        // Try permuting the order of the vertices
        list rStrat = permuteCompStrat(web, compStrat);
        compStrat = rStrat;
        kill rStrat;
        
        // If it's still no good, bail
        if( !compilationStratVerify(web, compStrat) )
        {
            print("[linkCompile] Unable to find a valid compilation strategy, failing.");
            return();
        }
    }
    kill web;
    
    // STEP 1. Go through each state of the link and compile the corresponding web. Then
    // add the data to returnListWebs. 
    // The compilation of a state 1101 is stored in web1101.
    int i;
    int numStateDiagrams = size(stateDiagrams);
    int globalZ2Degree = -1;
    
    for( i = 1; i <= numStateDiagrams; i++ )
    {
        // Our current state diagram:
        intvec v = stateDiagrams[i];
        
        dbprint(printlevel, "");
        dbprint(printlevel, "[linkCompile] Looking at state [" + string(v) + "] which is number [[ " + string(i) + "/" + string(numStateDiagrams) + " ]].");
        
        string listName = "compiled" + intvecToString(v); // compilation product
        string listName2 = "web" + intvecToString(v); // web
        
        string cmd = "int definedCheck = defined(" + listName + ") && defined(" + listName2 + ");";
        execute(cmd);
        
        if( definedCheck )
        {
            // We have already compiled this web. /////////////
            dbprint(printlevel, "[linkCompile] Already compiled, skipping.");
        }
        else
        {
            // We have not already compiled this web. ///////////
            dbprint(printlevel, "[linkCompile] Currently used memory: " + string(memory(0) div 1024) + "kb."); 
        
            // Build a web suitable for calling webCompile:
            list web = webFromLink(linkData, v, N);
            execute("list " + listName2 + " = web;");
            
            // Call webCompile:
            execute("list " + listName + " = webCompile(web, compStrat);");
            
            // Save the results to the savefile:
            dbprint(printlevel, "[linkCompile] Saving to disk.");

            execute("writeObjectToFile(fileName,\"" + listName + "\"," + listName + ");");
            execute("writeObjectToFile(fileName,\"" + listName2 + "\"," + listName2 + ");");
            
            kill web;
        }
        
        string cmd = "list Z = " + listName + ";";
        execute(cmd);
        
        // The return value of webCompile should look like "nonsquare", r1,r2, Q
        // and the last entry of Q should be ("nonsquare",r1,r2,R) with the last entry of R
        // a tuple (f0,f1,g0,g1,u,N) with f0,g0 splitting the even part of the very last
        // idempotent and similarly for f1,g1.
        if( typeof(Z[1]) != "string" || size(Z) != 4 || typeof(Z[4]) != "list" )
        {
            print("[linkCompile] Return of webCompile is not of the expected format, exiting.");
            return();
        }
        
        // Q = Z[4], Q[size(Q)] is ("nonsquare",r1,r2,R) so R is
        list R = Z[4][size(Z[4])][4];
        list lastEntry = R[size(R)];
        
        matrix f(0) = lastEntry[2];
        matrix f(1) = lastEntry[3];
        matrix g(0) = lastEntry[4];
        matrix g(1) = lastEntry[5];
        
        // First compute the dimension of the compilation in each degree:
        int j;
        for(j=0; j<=1; j++)
        {
            // The compilation of the current web is a Z/2-graded Q-vector space Y = Y0 + Y1.
            // The product f(i) * g(i) will either be a unitmatrix of size k if dim(Yi) = k > 0
            // or a 1x1 zero matrix if Yi = 0.
            int dimComp(j);
            if( f(j) * g(j) != matrix(0,1,1) && f(j) * g(j) != unitmat(nrows(f(j))) )
            {
                print("[linkCompile] Unexpected output from webCompile, exiting.");
                return();
            }
            
            // Y is zero in this degree.
            if( f(j) * g(j) == matrix(0,1,1) )
            {
                dimComp(j) = 0;
            }
            
            // Y is nonzero in this degree.
            if( f(j) * g(j) == unitmat(nrows(f(j))) )
            {
                dimComp(j) = nrows(f(j));
            }
            
            kill f(j), g(j);
        }
        
        // Decide on the degree in which we are nonzero and set the rank in that degree
        int rankk, compDegree;
        intvec grading;
        if( dimComp(0) == 0 && dimComp(1) != 0 )
        {
            // We are concentrated in degree 1.
            compDegree = 1;
            rankk = dimComp(1);
            grading = Z[3];
        }
        
        if( dimComp(0) != 0 && dimComp(1) == 0 )
        {
            compDegree = 0;
            rankk = dimComp(0);
            grading = Z[2];
        }
        
        if( dimComp(0) == 0 && dimComp(1) == 0 )
        {
            dbprint(printlevel, "[mfCompileLink] Web " + intvecToString(v) + " compiled to zero.");
            compDegree = -1;
            rankk = 0;
        }
        
        // All webs compilations should be living in one of the Z2 degrees.
        if( compDegree != -1 )
        {
            if( globalZ2Degree != -1 && globalZ2Degree != compDegree )
            {
                print("[linkCompile] The compilation of " + intvecToString(v) + " is in Z2 degree " + string(compDegree));
                print("[linkCompile] whereas an earlier web was in degree " + string(globalZ2Degree) + ", exiting.");
                return();
            }
            
            if( globalZ2Degree == -1 )
            {
                globalZ2Degree = compDegree;
            }
        }
        
        // Add the compilation products to the return list.
        // Recall that tuples in this list have the format
        // (intvec state, int rankk, int degree, intvec grading).
        list compiledEntry = v, rankk, compDegree, grading;
        returnListWebs = returnListWebs + list(compiledEntry);
        kill compiledEntry, v, Z, R, lastEntry, rankk, compDegree, grading;
        kill dimComp(0), dimComp(1);
    }
    
    dbprint(printlevel, "[linkCompile] Time taken to compile webs: " + string(timer-timeLink) + "ms."); 
    
    // STEP 2. Compile the maps
    dbprint(printlevel, "");
    dbprint(printlevel, "[linkCompile] Now compiling maps.");
    
    for( i = 1; i <= numStateDiagrams; i++ )
    {
        // Our current state diagram:
        intvec v = stateDiagrams[i];
        dbprint(printlevel, "[linkCompile] Looking at state [" + string(v) + "] which is number [[ " + string(i) + "/" + string(numStateDiagrams) + " ]].");
    
        // Look at each position and the possible maps given by chis.
        int j;
        for( j = 1; j <= size(v); j++ )
        {
            // If either the current position in v is a 0 and the crossing is positive, or
            // if the current position is 1 and the crossing is negative, then there is a map
            // originating in this state diagram with a chi in this position.
            intvec w = v;
            int somethingToDo = 0;
            int chitype;
        
            if(v[j] == 0 && vertexParity[j] == 1) // chitype is 0
            {
                w[j] = w[j]+1;
                somethingToDo = 1;
                chitype = 0;
            }
        
            if(v[j] == 1 && vertexParity[j] == -1) // chitype is 1
            {
                w[j] = w[j]-1;
                somethingToDo = 1;
                chitype = 1;
            }
        
            if( somethingToDo )
            {
                dbprint(printlevel, "[linkCompile]  Looking at the map to [" + string(w) + "].");
                
                // Load up the appropriate webs:
                string cmd = "list leftCompiledWeb = compiled" + intvecToString(v) + ";"; execute(cmd);
                cmd = "list leftWeb = web" + intvecToString(v) + ";"; execute(cmd);
                cmd = "list rightCompiledWeb = compiled" + intvecToString(w) + ";"; execute(cmd);
                cmd = "list rightWeb = web" + intvecToString(w) + ";"; execute(cmd);
                
                // Lookup the degree in which these compilation products are nonzero.
                // At position i in returnListWebs is various information about the
                // compilation of the "left" web.
                list bleft = returnListWebs[i];
                
                if( bleft[1] != v )
                {
                    print("[linkCompile] Unexpected data format, exiting.");
                    return();
                }
                
                int concentrateDegree = bleft[3];
                
                // Check whether the target web w compiles to zero (in the sense that
                // its differential is zero, not necessarily its rank):
                for( int aa=1; aa<=size(returnListWebs); aa++ )
                {
                    if( returnListWebs[aa][1] == w )
                    {
                        int wPosition = aa;
                        break;
                    }
                }
                
                list bright = returnListWebs[wPosition];
                int concentrateDegreeW = bright[3];
                
                if( concentrateDegree == -1 || concentrateDegreeW == -1 )
                {
                    // The source web compiles to zero (in the sense that its differential 
                    // is zero, not necessarily its rank), so we can just put some
                    // zero matrix onto the return list and skip.
                    
                    // In listWebsByDegree (and then in KRcomplex) we discard
                    // webs compiling to zero, so we should also not include zero maps
                    // from or to zero modules.
                    //
                    // list compiledEntry = v, w, v[j], j, zeromat(0,1,1);
                    // returnListMaps = returnListMaps + list(compiledEntry);
                    kill cmd, compiledEntry, leftCompiledWeb, rightCompiledWeb, leftWeb, rightWeb;

                    j++;
                    continue;
                }
                
                // Otherwise we have a well-defined degree, and a well-defined vector space
                // of the following rank in the source web.
                int rankSource = bleft[2];
                
                // Now check to see if we have already compiled this map:
                string listName = "map" + intvecToString(v) + "to" + intvecToString(w);

                string cmd = "int definedCheck = defined(" + listName + ");";
                execute(cmd);
                                
                if( definedCheck )
                {
                    // We have already compiled this morphism.
                    dbprint(printlevel, "[linkCompile]  Already compiled, skipping.");
                    
                    string cmd = "list D = " + listName + ";";
                    execute(cmd);
                }
                else
                {
                    // We have not already compiled this morphism.
                    intvec e = edgeList[j];
                    matrix M = mfChi(chitype, var(e[1]), var(e[2]), var(e[3]), var(e[4]), N);
                    kill e;
                
                    def Dp = webCompileMorphismFromData( leftWeb, leftCompiledWeb, rightWeb, rightCompiledWeb, compStrat, M, j );
                    
                    // We should end up with a zero potential, so the output of webCompileMorphism should be a pair,
                    // giving matrices in odd and even degrees
                    if( typeof(Dp) != "list" )
                    {
                        print("[linkCompile] webCompileMorphism did not return a pair, exiting.");
                        return();
                    }
                    list D = Dp;
                    kill Dp;
                    
                    string cmd = "list " + listName + " = D;"; execute(cmd);
                    kill M;
                    
                    // Check that D is of degree +1:
                    if( useSanityChecks )
                    {
                        matrix Dc = D[concentrateDegree+1];
                        int i2, j2;
                        for(i2 = 1; i2 <= nrows(Dc); i2++) // row
                        {
                            for(j2 = 1; j2 <= ncols(Dc); j2++) // column
                            {
                                if( Dc[i2,j2] != 0 )
                                {
                                    if( 1 - bright[4][i2] + bleft[4][j2] != 2 * deg(Dc[i2,j2]) )
                                    {
                                        print("[linkCompile] !!!WARNING!!! Map is not degree +1.");
                                        print("[linkCompile] In position (" + string(i2) + "," + string(j2) + ") there is " + string(Dc[i2,j2]) + "." +
                                                "The source has degree " + string(bleft[4][j2]) + " and the target " + string(bright[4][i2]) + ".");
                                    }
                                }
                            }
                        }
                        kill Dc;
                    }
                    
                    // Save the results to the savefile:
                    dbprint(printlevel, "[linkCompile] Saving to disk.");
                    execute("writeObjectToFile(fileName,\"" + listName + "\",D);");
                }
                
                if( ncols(D[concentrateDegree+1]) != rankSource )
                {
                    print("[linkCompile] Size is mismatch between map (" + string(nrows(D[concentrateDegree+1])) + "x" + string(ncols(D[concentrateDegree+1])) + ") and its source vector space (dim " + string(rankSource) + "), exiting.");
                    return();
                }
                
                // Add the computed list to the end of compiledEntry:
                list compiledEntry = v, w, v[j], j, D[concentrateDegree+1];                
                
                returnListMaps = returnListMaps + list(compiledEntry);
                kill listName, cmd, D, compiledEntry, leftWeb, rightWeb, leftCompiledWeb, rightCompiledWeb;
            }
            
            kill w, somethingToDo;
        }
        
        kill v;
    }
    
    setring RRR;
   
    dbprint(printlevel, "[linkCompile] total elapsed time " + string(timer-timeLink) + "ms."); 
    return(list(returnListWebs,returnListMaps));
}


////////////////////////////////////////////////////////////////////
// linkCompileParts
//
// Just like linkCompile, but computes only states (and not their morphisms) 
// in the range [startpos, endpos].

// TODO: Either remove or update to bug fixes in linkCompile and linkCompileReduced

proc linkCompileParts(list linkData, string fileName, int N, int startpost, int endpos)
{
    // We do not change rings: we assume that the ambient ring has the same number of variables
    // as edges in the link, and we associate ith variable var(i) with the edge i. Of course
    // if there is a save file, then loading this save file will set the ring.
    int useSanityChecks = blowFlags("sanity_checks"); 
    int useMPfile = blowFlags("use_mpfile");
    
    system("--ticks-per-sec",1000);
    int timeLink = timer;
 
    def RRR = basering;
    
    intvec vertexParity = linkData[1];
    list edgeList = linkData[2];
    int webSiz = size(vertexParity);
    
    // Read in the save file, if it exists
    if( status(fileName,"exists") == "yes" )
    {
        if( useMPfile )
        {
            def nR = loadSaveFile(fileName);
            setring nR;
        }
        else
        {
            execute(loadSaveFile(fileName));
        }
    }
        
    list returnListWebs;
    list returnListMaps;
    
    // The compilation strategy is the same for all webs: we simply compile from the 1st
    // vertex to the last one, and enumerate the variables in the given order (note that
    // all the ring variables are taken as internal variables).
    list compStratVert;
    list compStratVar;
    int j;
    for(j=1;j<=size(edgeList); j++)
    {
        compStratVert[j] = j;
    }
    for(j=1;j<=nvars(basering);j++)
    {
        compStratVar = compStratVar + list(var(j));
    }
    list compStrat = compStratVert,compStratVar;
        
    // STEP 1. Go through each state of the link and compile the corresponding web. Then
    // add the data to returnListWebs. 
    // The compilation of a state 1101 is stored in web1101.
    int i;
    list stateDiagrams = powervec(webSiz);
    int numStateDiagrams = size(stateDiagrams);
    
    int globalZ2Degree = -1;
    
    for( i = startpost; i <= endpos; i++ )
    {
        // Our current state diagram.
        intvec v = stateDiagrams[i];
        
        dbprint(printlevel, "");
        dbprint(printlevel, "[linkCompile] Looking at state [" + string(v) + "] which is number [[ " + string(i) + "/" + string(numStateDiagrams) + " ]].");
        
        string listName = "compiled" + intvecToString(v); // compilation product
        string listName2 = "web" + intvecToString(v); // web
        
        string cmd = "int definedCheck = defined(" + listName + ") && defined(" + listName2 + ");";
        execute(cmd);
        
        if( definedCheck )
        {
            // We have already compiled this web. /////////////
            dbprint(printlevel, "[linkCompile] Already compiled, skipping.");
        }
        else
        {
            // We have not already compiled this web. ///////////
            dbprint(printlevel, "[linkCompile] Currently used memory: " + string(memory(0) div 1024) + "kb."); 
        
            // Build a web suitable for calling webCompile:
            list web = webFromLink(linkData, v, N);
            
            execute("list " + listName2 + " = web;");
            
            // Call webCompile:
            execute("list " + listName + " = webCompile(web, compStrat);");
            
            // Save the results to the savefile:
            dbprint(printlevel, "[linkCompile] Saving to disk.");

            execute("writeObjectToFile(fileName,\"" + listName + "\"," + listName + ");");
            execute("writeObjectToFile(fileName,\"" + listName2 + "\"," + listName2 + ");");
            
            kill web;
        }
        
        string cmd = "list Z = " + listName + ";";
        execute(cmd);
        
        // The return value of webCompile should look like "nonsquare", r1,r2, Q
        // and the last entry of Q should be ("nonsquare",r1,r2,R) with the last entry of R
        // a tuple (f0,f1,g0,g1,u,N) with f0,g0 splitting the even part of the very last
        // idempotent and similarly for f1,g1.
        if( typeof(Z[1]) != "string" || size(Z) != 4 || typeof(Z[4]) != "list" )
        {
            print("[linkCompile] Return of webCompile is not of the expected format, exiting.");
            return();
        }
        
        // Q = Z[4], Q[size(Q)] is ("nonsquare",r1,r2,R) so R is
        list R = Z[4][size(Z[4])][4];
        list lastEntry = R[size(R)];
        
        matrix f(0) = lastEntry[2];
        matrix f(1) = lastEntry[3];
        matrix g(0) = lastEntry[4];
        matrix g(1) = lastEntry[5];
        
        // First compute the dimension of the compilation in each degree.
        int j;
        for(j=0; j<=1; j++)
        {
            // The compilation of the current web is a Z/2-graded Q-vector space Y = Y0 + Y1.
            // The product f(i) * g(i) will either be a unitmatrix of size k if dim(Yi) = k > 0
            // or a 1x1 zero matrix if Yi = 0.
            int dimComp(j);
            if( f(j) * g(j) != matrix(0,1,1) && f(j) * g(j) != unitmat(nrows(f(j))) )
            {
                print("[linkCompile] Unexpected output from webCompile, exiting.");
                return();
            }
            
            // Y is zero in this degree
            if( f(j) * g(j) == matrix(0,1,1) )
            {
                dimComp(j) = 0;
            }
            
            // Y is nonzero in this degree
            if( f(j) * g(j) == unitmat(nrows(f(j))) )
            {
                dimComp(j) = nrows(f(j));
            }
            
            kill f(j), g(j);
        }
        
        // Decide on the degree in which we are nonzero and set the rank in that degree.
        int rankk, compDegree;
        intvec grading;
        if( dimComp(0) == 0 && dimComp(1) != 0 )
        {
            // We are concentrated in degree 1
            compDegree = 1;
            rankk = dimComp(1);
            grading = Z[3];
        }
        
        if( dimComp(0) != 0 && dimComp(1) == 0 )
        {
            compDegree = 0;
            rankk = dimComp(0);
            grading = Z[2];
        }
        
        if( dimComp(0) == 0 && dimComp(1) == 0 )
        {
            dbprint(printlevel, "[mfCompileLink] Web " + intvecToString(v) + " compiled to zero.");
            compDegree = -1;
            rankk = 0;
        }
        
        // All webs compilations should be living in one of the Z2 degrees.
        if( compDegree != -1 )
        {
            if( globalZ2Degree != -1 && globalZ2Degree != compDegree )
            {
                print("[linkCompile] The compilation of " + intvecToString(v) + " is in Z2 degree " + string(compDegree));
                print("[linkCompile] whereas an earlier web was in degree " + string(globalZ2Degree) + ", exiting.");
                return();
            }
            
            if( globalZ2Degree == -1 )
            {
                globalZ2Degree = compDegree;
            }
        }
        
        // Add the compilation products to the return list.
        // Recall that tuples in this list have the format
        // (intvec state, int rankk, int degree, intvec grading)
        list compiledEntry = v, rankk, compDegree, grading;
        returnListWebs = returnListWebs + list(compiledEntry);
        kill compiledEntry, v, Z, R, lastEntry, rankk, compDegree, grading;
        kill dimComp(0), dimComp(1);
    }
    
    dbprint(printlevel, "[linkCompile] Time taken to compile webs: " + string(timer-timeLink) + "ms."); 

    setring RRR;
   
    dbprint(printlevel, "[linkCompile] total elapsed time " + string(timer-timeLink) + "ms."); 
    kill returnListWebs;
    return(list());
}




////////////////////////////////////////////////////////////////////
// linkCompileReduced
//
// Almost identical to linkCompile, except that we take an additional
// param (a ring variable reductionMark) and use it to compute the reduced
// bigraded sl(N) Khovanov-Rozansky homology. Essentially we take the webs from
// the usual calculation, drop the edge with reductionMark on it, and substitute
// reductionMark = 0 in all MFs.
//
// The return is a pair (X,Y) of lists:
//
//   1. X is a list of tuples (intvec state, int rank, -1, intvec grading) where state ranges
//   over all state diagrams, and grading is the grading vector, rank = size(grading).
//
//   2. Y is a list of tuples (intvec state1,intvec state2,int chitype,int chipos,matrix F).
//   where state1, state2 range over all state diagrams, chitype is 0 or 1 resp. if the map
//   is a chi0 or chi1 map, and chipos gives the position in the state diagram where the
//   map is inserted. Finally, F is the matrix of the map, the result of calling
//   webCompileMorphism.
//
// NOTE: The final param tells us which part of the Z/2-grading to use (0 for even, 1 for
// odd, and -1 for both).

proc linkCompileReduced(list linkData, string fileName, int N, poly reductionMark, int useDegree)
{
    // We do not change rings: we assume that the ambient ring has the same number of variables
    // as edges in the link, and we associate ith variable var(i) with the edge i. Of course
    // if there is a save file, then loading this save file will set the ring.
    int useSanityChecks = blowFlags("sanity_checks"); 
    int useMPfile = blowFlags("use_mpfile");
    
    if( useDegree != 0 && useDegree != 1 && useDegree != -1 )
    {
        print("[linkCompileReduced] useDegree should be in {0,1,-1}, exiting.");
        return();
    }
    
    system("--ticks-per-sec",1000);
    int timeLink = timer;
 
    def RRR = basering;
    
    intvec vertexParity = linkData[1];
    list edgeList = linkData[2];
    int webSiz = size(vertexParity);
    
    // Read in the save file, if it exists:
    if( status(fileName,"exists") == "yes" )
    {
        if( useMPfile )
        {
            def nR = loadSaveFile(fileName);
            setring nR;
            poly reductionMark = imap(RRR,reductionMark);
        }
        else
        {
            execute(loadSaveFile(fileName));
        }
    }
        
    list returnListWebs;
    list returnListMaps;
    
    // The compilation strategy is the same for all webs. We always enumerate the variables
    // in the default order, but we have to choose an ordering of the vertices which gives
    // a valid compilation strategy. That is, after fusing the first i vertices, the result
    // should be connected to the (i+1)st vertex by an edge. Note that we delete the edge
    // with the variable reductionMark, so this edge doesn't count!
    
    // First build the default compilation strategy
    list compStratVert;
    list compStratVar;
    int j;
    for(j=1;j<=size(edgeList); j++)
    {
        compStratVert[j] = j;
    }
    for(j=1;j<=nvars(basering);j++)
    {
        // Do not add reductionMark as an internal variable
        if( reductionMark != var(j) )
        {
            compStratVar = compStratVar + list(var(j));
        }
    }
    list compStrat = compStratVert,compStratVar;
        
    // Now make sure it is valid, and permute the ordering of the vertices if not
    list stateDiagrams = powervec(webSiz);
    list web = webFromLink(linkData, stateDiagrams[1], N, reductionMark);
    
    if( !compilationStratVerify(web, compStrat) )
    {
        dbprint(printlevel, "[linkCompileReduced] Default compilation strategy is invalid, trying to permute vertices...");
        
        // Try permuting the order of the vertices
        list rStrat = permuteCompStrat(web, compStrat);
        compStrat = rStrat;
        kill rStrat;
        
        // If it's still no good, bail
        if( !compilationStratVerify(web, compStrat) )
        {
            print("[linkCompileReduced] Unable to find a valid compilation strategy, failing.");
            return();
        }
    }
    kill web;
    
    // STEP 1. Go through each state of the link and compile the corresponding web. Then
    // add the data to returnListWebs. 
    // The compilation of a state 1101 is stored in web1101.
    int i;
    int numStateDiagrams = size(stateDiagrams);
    
    int globalZ2Degree = -1;
    
    for( i = 1; i <= numStateDiagrams; i++ )
    {
        // Our current state diagram:
        intvec v = stateDiagrams[i];
        
        dbprint(printlevel, "");
        dbprint(printlevel, "[linkCompileReduced] Looking at state [" + string(v) + "] which is number [[ " + string(i) + "/" + string(numStateDiagrams) + " ]].");
        
        string listName = "compiled" + intvecToString(v); // compilation product
        string listName2 = "web" + intvecToString(v); // web
        
        string cmd = "int definedCheck = defined(" + listName + ") && defined(" + listName2 + ");";
        execute(cmd);
        
        if( definedCheck )
        {
            // We have already compiled this web. /////////////
            dbprint(printlevel, "[linkCompileReduced] Already compiled, skipping.");
        }
        else
        {
            // We have not already compiled this web. ///////////
            dbprint(printlevel, "[linkCompileReduced] Currently used memory: " + string(memory(0) div 1024) + "kb."); 
        
            // Build a web suitable for calling webCompile:
            list web = webFromLink(linkData, v, N, reductionMark);
            execute("list " + listName2 + " = web;");
            
            // Call webCompile:
            execute("list " + listName + " = webCompile(web, compStrat);");
            
            // Save the results to the savefile:
            dbprint(printlevel, "[linkCompileReduced] Saving to disk.");

            execute("writeObjectToFile(fileName,\"" + listName + "\"," + listName + ");");
            execute("writeObjectToFile(fileName,\"" + listName2 + "\"," + listName2 + ");");
            
            kill web;
        }
        
        string cmd = "list Z = " + listName + ";";
        execute(cmd);
        
        // In contrast to the nonreduced homology, the return value of webCompile may live
        // in more than one Z2 degree. In fact we assume that it does.
        
        if( typeof(Z[1]) != "matrix" || typeof(Z[2]) != "intvec" || size(Z[2]) != ncols(Z[1]) )
        {
            print("[linkCompileReduced] Return of webCompile is not of the expected format, exiting.");
            return();
        }
        
        // We also assume that Z[1] is a zero matrix
        if( Z[1] != matrix(0,nrows(Z[1]),ncols(Z[1])) )
        {
            print("[linkCompileReduced] Return of webCompile not a zero matrix, exiting.");
            return();
        }
        
        // Add the compilation products to the return list
        // Recall that tuples in this list have the format
        // (intvec state, int rank, int concdegree, intvec grading)
        // We match the format of linkCompile just for convenience, although of
        // course the concentration degree isn't relevant here.
        
        if( useDegree == -1 )
        {
            list compiledEntry = v, size(Z[2]), -1, Z[2];
        }
        
        if( useDegree == 1 )
        {
            intvec b = Z[2];
            intvec Zodd = b[(size(Z[2]) div 2 + 1)..size(Z[2])];
            list compiledEntry = v, size(Zodd), -1, Zodd;
            kill b, Zodd;
        }
        
        if( useDegree == 0 )
        {
            intvec b = Z[2];
            intvec Zeven = b[1..(size(Z[2]) div 2)];
            list compiledEntry = v, size(Zeven), -1, Zeven;
            kill b, Zeven;
        }
        
        returnListWebs = returnListWebs + list(compiledEntry);
        kill compiledEntry, v, Z;
    }
    
    dbprint(printlevel, "[linkCompileReduced] Time taken to compile webs: " + string(timer-timeLink) + "ms."); 
    
    // STEP 2. Compile the maps.
    dbprint(printlevel, "");
    dbprint(printlevel, "[linkCompileReduced] Now compiling maps.");
    
    for( i = 1; i <= numStateDiagrams; i++ )
    {
        // Our current state diagram:
        intvec v = stateDiagrams[i];
        dbprint(printlevel, "[linkCompileReduced] Looking at state [" + string(v) + "] which is number [[ " + string(i) + "/" + string(numStateDiagrams) + " ]].");
    
        // Look at each position and the possible maps given by chis.
        int j;
        for( j = 1; j <= size(v); j++ )
        {
            // If either the current position in v is a 0 and the crossing is positive, or
            // if the current position is 1 and the crossing is negative, then there is a map
            // originating in this state diagram with a chi in this position.
            intvec w = v;
            int somethingToDo = 0;
            int chitype;
        
            if(v[j] == 0 && vertexParity[j] == 1) // chitype is 0
            {
                w[j] = w[j]+1;
                somethingToDo = 1;
                chitype = 0;
            }
        
            if(v[j] == 1 && vertexParity[j] == -1) // chitype is 1
            {
                w[j] = w[j]-1;
                somethingToDo = 1;
                chitype = 1;
            }
        
            if( somethingToDo )
            {
                dbprint(printlevel, "[linkCompileReduced]  Looking at the map to [" + string(w) + "].");
                
                // Load up the appropriate webs:
                string cmd = "list leftCompiledWeb = compiled" + intvecToString(v) + ";"; execute(cmd);
                cmd = "list leftWeb = web" + intvecToString(v) + ";"; execute(cmd);
                cmd = "list rightCompiledWeb = compiled" + intvecToString(w) + ";"; execute(cmd);
                cmd = "list rightWeb = web" + intvecToString(w) + ";"; execute(cmd);
                
                // Now check to see if we have already compiled this map. Note that
                // in distinction to the unreduced case, the output of webCompileMorphism
                // is a matrix not a list!
                string matrixName = "map" + intvecToString(v) + "to" + intvecToString(w);

                string cmd = "int definedCheck = defined(" + matrixName + ");";
                execute(cmd);
                                
                if( definedCheck )
                {
                    // We have already compiled this morphism.
                    dbprint(printlevel, "[linkCompileReduced]  Already compiled, skipping.");
                    
                    string cmd = "matrix D = " + matrixName + ";";
                    execute(cmd);
                }
                else
                {
                    // We have not already compiled this morphism.
                    intvec e = edgeList[j];
                    matrix M = mfChi(chitype, var(e[1]), var(e[2]), var(e[3]), var(e[4]), N);
                    M = subst(M,reductionMark,0);
                    kill e;

                    def Dp = webCompileMorphismFromData( leftWeb, leftCompiledWeb, rightWeb, rightCompiledWeb, compStrat, M, j );
                    
                    // The output of webCompileMorphism should be a matrix.
                    if( typeof(Dp) != "matrix" )
                    {
                        print("[linkCompileReduced] webCompileMorphism did not return a matrix.");
                        return();
                    }
                    matrix D = Dp;
                    kill Dp;
                    
                    // Other sanity checks:
                    if( ncols(D) != ncols(leftCompiledWeb[1]) || nrows(D) != ncols(rightCompiledWeb[1]) )
                    {
                        print("[linkCompileReduced] Size is mismatch, exiting.");
                        return();
                    }
                    
                    string cmd = "matrix " + matrixName + " = D;"; execute(cmd);
                    kill M;

                    // Save the results to the savefile:
                    dbprint(printlevel, "[linkCompileReduced] Saving to disk.");
                    execute("writeObjectToFile(fileName,\"" + matrixName + "\",D);");
                }
                
                if( useDegree == -1 )
                {
                    list compiledEntry = v, w, v[j], j, D;
                }
                
                if( useDegree == 1 )
                {
                    list compiledEntry = v, w, v[j], j, extractblockmat(D)[4];
                }
                
                if( useDegree == 0 )
                {
                    list compiledEntry = v, w, v[j], j, extractblockmat(D)[1];
                }
                
                returnListMaps = returnListMaps + list(compiledEntry);
                kill matrixName, cmd, D, compiledEntry, leftWeb, rightWeb, leftCompiledWeb, rightCompiledWeb;
            }
            
            kill w, somethingToDo;
        }
        
        kill v;
    }
    
    setring RRR;
   
    dbprint(printlevel, "[linkCompileReduced] total elapsed time " + string(timer-timeLink) + "ms."); 
    return(list(returnListWebs,returnListMaps));
}

////////////////////////////////////////////////////////////////////
// loadSaveFile
//
// Given the filename of an MP or ASCII data file, reads in the contents. We require
// that all writes to the MP file be made in the context of a particular ring R.
// Then loadFromMPfile will create its own copy of this ring, and the loaded data
// types will belong to this ring. This ring is the single return value, in the case
// where MPfile is being used.
//
// So in this case loadSaveFile should be called as follows:
//
//   def nR = loadFromMPfile("example.mp");
//   setring nR;
//
// If we are in ASCII mode then we simply return the contents of the given
// file, so that running execute should recreate the needed objects in memory, and we
// use the current basering (there is no return value in this case).
//
// In either case objects should be written to the store using writeObjectToFile.

proc loadSaveFile(string filename)
{
    int useMPfile = blowFlags("use_mpfile");

    if( !status(filename,"exists","yes") )
    {
        print("[loadSaveFile] File does not exist, exiting.");
        return();
    }
    
    if( !useMPfile )
    {
        // ASCII mode
        return(read(filename));
    }
    
    // Now we are in MP mode. First open the link:
    link l = "MPfile:r " + filename;
    open(l);
    
    if( !status(l,"open","yes") )
    {
        print("[loadFromMPfile] Failed to open link, exiting.");
        return();
    }
    
    while( status(l,"read","ready") )
    {
        // Read the next line:
        def L = read(l);
        
        if( typeof(L) == "list" )
        {
            // The first position is the variable name, the second the actual object
            string cmd = "def " + L[1] + " = L[2];";
            execute(cmd);

            cmd = "exportto(Top," + L[1] + ");";
            execute(cmd);
            kill cmd;
        }
        
        kill L;
     }
     
     close(l);
     return(basering);
}
    
////////////////////////////////////////////////////////////////////
// writeObjectToFile
//
// Given an object and its name, we write it out to the given file.

proc writeObjectToFile(string fileName, string nam, def L)
{
    int useMPfile = blowFlags("use_mpfile");
    
    if( useMPfile )
    {
        // If we are using MPfile then we write out the list (nam, L).
        link l = "MPfile: " + fileName;
        execute("write(l,list(\"" + nam + "\",L));");
        close(l);
    }
    else
    {
        // If we are using ASCII, write out to the file the lines
        // list nam = ???; (or matrix nam = ??, etc) 
        // where ?? is the output of archiveListToString.
        string s = typeof(L) + " " + nam + " = " + archiveObjectToString(L) + ";";
        write(":a " + fileName, s);
    }
}

////////////////////////////////////////////////////////////////////
// archiveObjectToString
//
// Given an object, creates a string representation which allows us to
// reconstruct the list in memory exactly as it was before. There is
// a limitation to the kinds of lists that we allow.
//
// Kudos: this is inspired by writelist, part of inout.lib in the
// Singular standard libraries.

proc archiveObjectToString(def L)
{
    string ent;
    int foundType;
    
    if( typeof(L) == "int" )
    {
        ent = "int(" + string(L) + ")";
        foundType = 1;
    }
    
    if( typeof(L) == "number" )
    {
        ent = "number(" + string(L) + ")";
        foundType = 1;
    }
    
    if( typeof(L) == "string" )
    {
        ent = "string(\"" + L + "\")";
        foundType = 1;
    }
    
    if( typeof(L) == "intvec" )
    {
        ent = "intvec(" + string(L) + ")";
        foundType = 1;
    }
        
    if( typeof(L) == "matrix" )
    {
        ent = "matrix(ideal(" + string(L) + ")," + string(nrows(L)) + "," + string(ncols(L)) + ")";
        foundType = 1;
    }
        
    if( typeof(L) == "poly" )
    {
        ent = "poly(" + string(L) + ")";
        foundType = 1;
    }
        
    if( typeof(L) == "list" )
    {
        ent = "list(";
        int j;
        for(j=1;j<=size(L);j++)
        {
            if( j != 1 ){ ent = ent + ","; }
                
            ent = ent + archiveObjectToString(L[j]);
        }
            
        ent = ent + ")";
        foundType = 1;
    }
        
    if( !foundType )
    {
        print("[archiveObjectToString] Found unsupported type \"" + typeof(L) + "\" in list, skipping!");
    }
    
    return(ent);
}

/////////////////////////////////////////////////
// DATA TYPE: complex
//
// Let X be a complex of free modules (possibly of zero rank). We represent this
// as a pair (m, R, L[,G]) consisting of an integer m, a list R of non-negative integers
// and a list L of matrices. We should have size(R) = size(L) + 1.
//
// Here the rank of X^i should be zero for i < m. For i >= m either X^i is rank zero
// or it has positive rank. R[1] should be the rank of X^m, R[2] the rank of X^(m+1), and so on.
// L should be a list of matrices giving the differentials in X, L[1] the map X^m --> X^(m+1)
// and so on. In case a differential ends or begins on a free module of rank zero, the entry
// in L should be a 1x1 zero matrix.
//
// Optionally a third list G can be given, with G[1] an intvec giving the grading vector
// on the complex in degree m, and so on. If X^i is rank zero, then the corresponding entry
// of G is just intvec(0).
//
// Complexes should be accessed using the setter/getter methods below.
//
// EXAMPLE:
//
// The list (0,(1,1,1),(x,y)) describes the cpx R -- x --> R -- y --> R (say over R = k[x,y]/xy!)
// beginning in degree zero, whereas (0,(1,0,1),(0,0)) is R ---> 0 ---> R.

proc cpxMinDegree(list cpx){ return(cpx[1]); }

// Degree where the last differential ends:
proc cpxMaxDegree(list cpx){ return(cpx[1] + size(cpx[3])); }

proc cpxDifferentialInDegree(list cpx, int degr)
{
    int mind = cpxMinDegree(cpx);
    int maxd = cpxMaxDegree(cpx);
    
    if( degr < mind || degr > maxd - 1 )
    {
        print("[cpxDifferentialInDegree] Degree outside bounds, exiting.");
        return();
    }
    
    return(cpx[3][degr-mind+1]);
}

proc cpxRankInDegree(list cpx, int degr)
{
    // We deduce the rank from the number of rols/cols in the differential.
    int mind = cpxMinDegree(cpx);
    int maxd = cpxMaxDegree(cpx);
    
    if( degr < mind || degr > maxd ){ return(0); }

    return(cpx[2][degr-mind+1]);
}

proc cpxGradingInDegree(list cpx, int degr)
{
    if( !cpxIsGraded(cpx) )
    {
        print("[cpxGradingInDegree] Complex is not graded, exiting.");
        return();
    }
    
    int mind = cpxMinDegree(cpx);
    int maxd = cpxMaxDegree(cpx);
    
    if( degr < mind || degr > maxd )
    {
        print("[cpxGradingInDegree] Degree is outside bounds, exiting.");
        return();
    }
    
    return(cpx[4][degr-mind+1]);
}

proc cpxIsGraded(list cpx){ return (size(cpx) > 3); }

////////////////////////////////////////////////////////////////////
// cpxIsValid
//

proc cpxIsValid(list cpx)
{
    if( size(cpx[2]) != size(cpx[3]) + 1)
    {
        print("[cpxIsValid] Relative size of rank and differential lists is wrong, failed.");
        return(0);
    }
    
    int mind = cpxMinDegree(cpx);
    int maxd = cpxMaxDegree(cpx);
    
    int i;
    for(i=mind; i<maxd; i++)
    {
        matrix D = cpxDifferentialInDegree(cpx,i);
        if( ncols(D) != cpxRankInDegree(cpx,i) && cpxRankInDegree(cpx,i) != 0 )
        {
            print("[cpxIsValid] Wrong column number in differential " + string(i));
            return(0);
        }
        
        if( nrows(D) != cpxRankInDegree(cpx,i+1) && cpxRankInDegree(cpx,i+1) != 0 )
        {
            print("[cpxIsValid] Wrong row number in differential " + string(i));
            return(0);
        }
    }
    
    int i;
    for(i=mind+1; i<maxd; i++)
    {
        matrix D = cpxDifferentialInDegree(cpx,i);
        matrix E = cpxDifferentialInDegree(cpx,i-1);
        
        if( D * E != matrix(0,nrows(D),ncols(E)) )
        {
            print("[cpxIsValid] The composition d[" + string(i) + "] d[" + string(i-1) + "] is nonzero, failed.");
            return(0);
        }
    }
    
    // If grading is on, check the number of grades matches the rank
    if( cpxIsGraded(cpx) )
    {
        for(i=mind; i<=maxd; i++)
        {
            if( cpxRankInDegree(cpx,i) != 0 && size(cpxGradingInDegree(cpx,i)) != cpxRankInDegree(cpx,i) )
            {
                print("[cpxIsValid] Rank does not match size of grading vector, failed.");
                return(0);
            }
        }
    }
    
    return(1);
}

////////////////////////////////////////////////////////////////////
// complexTensor
//
// Given two complexes (possibly graded) return the tensor product.
//
// NOTE: If we are given gradings, we do not check that the differentials
// are degree zero.

proc cpxTensor(list cpx(1), list cpx(2))
{
    if( !cpxIsValid(cpx(1)) || !cpxIsValid(cpx(2)) )
    {
        print("[cpxTensor] Passed invalid complexes, exiting.");
        return();
    }

    int useGrading;
    if( cpxIsGraded(cpx(1)) ){ int useGrading = 1; }
    
    int i,j;
    for(i=1; i<=2; i++)
    {
        int minDegree(i) = cpxMinDegree(cpx(i));
        int maxDegree(i) = cpxMaxDegree(cpx(i));
    }
    
    int minDegreeTensor = minDegree(1) + minDegree(2);
    int maxDegreeTensor = maxDegree(1) + maxDegree(2);
    
    // Decide on the ranks of the free modules in the tensor product, and at
    // the same time record which tensors contribute to each degree of the product.
    list tensorRanks;
    list nonzeroTensors;
    if( useGrading ){ list nonzeroTensorsGrading; }
    
    for(i=1;i<=maxDegreeTensor-minDegreeTensor+1;i++)
    {
        tensorRanks[i] = 0;
        nonzeroTensors[i] = list();
        if( useGrading ){ nonzeroTensorsGrading[i] = list(); }
    }
    
    for(i=minDegree(1);i<=maxDegree(1);i++)
    {
        for(j=minDegree(2);j<=maxDegree(2);j++)
        {
            int delta = cpxRankInDegree(cpx(1),i) * cpxRankInDegree(cpx(2),j);
            
            if( delta != 0 )
            {
                nonzeroTensors[i+j-minDegreeTensor+1] = nonzeroTensors[i+j-minDegreeTensor+1] + list(intvec(i,j));
                
                if( useGrading )
                {
                    intvec gr(1) = cpxGradingInDegree(cpx(1),i);
                    intvec gr(2) = cpxGradingInDegree(cpx(2),j);
                    
                    intvec tensorgr;
                    int k,l;
                    int count;
                    for(k=1;k<=size(gr(1));k++)
                    {
                        for(l=1;l<=size(gr(2));l++)
                        {
                            count++;
                            tensorgr[count] = gr(1)[k] + gr(2)[l];
                        }
                    }
                    nonzeroTensorsGrading[i+j-minDegreeTensor+1] = nonzeroTensorsGrading[i+j-minDegreeTensor+1] + list(tensorgr);
                }
            }
            
            tensorRanks[i+j-minDegreeTensor+1] = tensorRanks[i+j-minDegreeTensor+1] + delta;
        }
    }
    
    // Next construct the differentials. Note that in the matrices we ignore free modules
    // of rank zero. Suppose we have free modules X of rank r and Y of rank s, and index
    // the summands of X as X1,...,Xr and of Y as Y1,...,Ys. Then the tensor product X x Y
    // is a direct sum of Xi x Yj, and we order these according to the lexicographic ordering
    // on pairs (i,j). I.e. (X1 + X2) x (Y1 + Y2) = X1 x Y1 + X1 x Y2 + X2 x Y1 + X2 x Y2.
    int k;
    list finalDiffs;
    for(k=minDegreeTensor; k<maxDegreeTensor; k++)
    {
        matrix D;
        int sk = k - minDegreeTensor + 1;
        
        //print("Constructing differential in degree " + string(k));
        
        // We are trying to construct the differential in degree k of the tensor product.
        if( size(nonzeroTensors[sk]) == 0 || size(nonzeroTensors[sk+1]) == 0 )
        {
            // The source or target has rank zero
            D = matrix(0,1,1);
        }
        else
        {
            list maps;
            for(j=1; j<=size(nonzeroTensors[sk+1]); j++)
            {
                for(i=1; i<=size(nonzeroTensors[sk]); i++)
                {
                    intvec a = nonzeroTensors[sk][i];
                    intvec b = nonzeroTensors[sk+1][j];
                 
                    //print("Looking at map from " + string(a) + " to " + string(b));
                    
                    // say a = (a1,a2) and b = (b1,b2). Then the "a summand" of the tensor
                    // product in degree a1 + a2 = k has rank
                    int rk_src = cpxRankInDegree(cpx(1),a[1]) * cpxRankInDegree(cpx(2),a[2]);
                    int rk_tar = cpxRankInDegree(cpx(1),b[1]) * cpxRankInDegree(cpx(2),b[2]);
                    
                    matrix Z[rk_tar][rk_src];
                    
                    if( b[1] == a[1] + 1 && a[2] == b[2] )
                    {
                        // The map (a[1],a[2]) --> (a[1] + 1, a[2])
                        Z = tensor(cpxDifferentialInDegree(cpx(1),a[1]), unitmat(cpxRankInDegree(cpx(2),b[2])));
                    }
                    if( b[1] == a[1] && b[2] == a[2] + 1 )
                    {
                        // The map (a[1],a[2]) --> (a[1], a[2]+1)
                        Z = (-1)^(a[1]) * tensor(unitmat(cpxRankInDegree(cpx(1),a[1])), cpxDifferentialInDegree(cpx(2),a[2]));
                    }
                    
                    //print(Z);
                    //print("----");
                    maps = maps + list(Z);
                }
            }
            
            // Assemble all the maps
            D = superblockmat(size(nonzeroTensors[sk]),maps);
        }
        
        finalDiffs[sk] = D;
    }
    
    if( useGrading )
    {
        // Construct the grading vectors:
        list finalGr;
        for(k=minDegreeTensor; k<=maxDegreeTensor; k++)
        {
            intvec grade;
            int sk = k - minDegreeTensor + 1;
            
            if( size(nonzeroTensorsGrading[sk]) > 0 )
            {
                grade = nonzeroTensorsGrading[sk][1];
                
                for(j=2; j<=size(nonzeroTensorsGrading[sk]); j++)
                {
                    grade = grade, nonzeroTensorsGrading[sk][j];
                }
            }
            finalGr = finalGr + list(grade);
            kill grade;
        }
    }
    
    list tensorCpx = minDegreeTensor, tensorRanks, finalDiffs;
    
    if( useGrading )
    {
        tensorCpx = tensorCpx + list(finalGr);
    }
    
    if( !cpxIsValid(tensorCpx) )
    {
        print("[cpxTensor] Failed to construct valid complex, exiting.");
        return();
    }
    
    return(tensorCpx);
}