zapgraph2 := proc(G)
# Puts selected graphs in normal form, and returns empty set
# if the graph is unsuitable for putting in normal form 
local i, j, k, currcount, currnum, lastcount, lastnum, L, H, searching, P;
   L := [[],[]];
   searching := true;
   i := 0;
   lastnum := 2;
   lastcount := 4;
   currcount := 0;
   currnum := 3;
   while searching do
     i := i + 1;
     if nops(G) >= i and G[i][1]=1 then
        if not(G[i][2]=2) then
           if G[i][2]=currnum then
              currcount := currcount + 1
           else
              if (G[i][2] > (currnum + 1)) or 
                   (currcount > lastcount) or (currcount > 2) then 
                     RETURN({}) else
                 L := subsop(currcount = [op(L[currcount]),currnum],L);
                 lastnum := currnum;
                 lastcount := currcount;
                 currcount := 1;
                 currnum := G[i][2]
              fi  
           fi
        fi
     else
        searching := false;
        if (currcount > lastcount) or (currcount > 2) then RETURN({}) else
           L := subsop(currcount = [op(L[currcount]), currnum], L) 
        fi
     fi
   od;
   L := [[op({op(vertices(G))} minus ({1,2} union {op(L[1])} 
        union {op(L[2])}))], op(L)];
   H := G;
   if not(L[3]=[]) then
      P := combinat[permute](L[3]);
      P := map(M -> [1,2,op(M),op(L[2]),op(L[1])], P);
      H := pnf(H,P)
   fi;
   if not(L[2]=[]) then
      P := combinat[permute](L[2]);
      P := map(M -> [1,2,op(L[3]), op(M),op(L[1])], P);
      H := pnf(H,P)
   fi;
   if not(L[1]=[]) then
      P := combinat[permute](L[1]);
      P := map(M -> [1,2,op(L[3]),op(L[2]),op(M)], P);
      H := pnf(H,P)
   fi;
   {H}
end;
   

notriangleq := proc(L)
# Check for a triangle in the graph L[1] not at the basepoint
# where the subforest L[2] contains exactly one edge of this 
# triangle.
# See Lemma 4.3 of [H-V]
local n, i, j, k, G, F;
   G := L[1];
    F := L[2];
    n := v(G);
    for i from 2 to n do for j from i + 1 to n do for k from j + 1 to n do
                if nops({[i, j], [i, k], [j, k]} intersect {op(G)}) = 3
                 and nops({[i, j], [i, k], [j, k]} intersect {op(F)}) = 1
                then RETURN(false)
                fi
            od
        od
    od;
    true
end;


nodoubleq := proc(G)
# Search for double edge not at basepoint in graph G. 
# See Lemma 4.2 of [H-V].
local H;
   H := G;
   while not(H=[]) and H[1][1] = 1 do
     H := H[2..nops(H)]
   od;
   if nops(H) = nops({op(H)}) then true else false fi 
end;

collapse := proc(G, M)
# Collapse edge M from graph G
local i, j, c, d, L, exc, unzapped, k;
    i := M[1];
    j := M[2];
    unzapped := true;
    k := 1;
    L := G;
    while unzapped and k <= nops(L) do
        if L[k][1] = i and L[k][2] = j then
            unzapped := false; L := subsop(k = NULL, L)
        fi;
        k := k + 1
    od;
    c := proc(k, i, j)
            if k < j then k
            elif k = j then i
            elif j < k and 2 < j then k - 1
            elif j < k and j = 2 then k
            fi
        end;
    L := select(M -> not(M[2]=1), 
          insort(map((M, i, j, c) -> map(c, M, i, j), L, i, j, c)))
end;

findplusfaces := proc(G, forests)
# Given the graph and the collection of subforests which give the
# plus faces, this returns a list of pairs, where the first part
# of each pair is the graph and the second is the relevant subforest
local plusfaces, i;
    plusfaces := [];
    for i to nops(forests) do
        plusfaces := [op(plusfaces), [G, forests[i]]]
    od;
    plusfaces
end;

findfandb := proc(nullplus, fminus, bm)
# Returns information in a quadrouple about how to take chains of
# cubes which map to zero under the boundary map only on the plusfaces.
# First part of the quadrouple is a list of coefficients for the cubes 
# in the chains.  Second part is a list of coefficients for the boundary
# map of the chains on the minusfaces.  The third and fourth parts,
# respectively, are the actual cube and minusface lists that correspond
# to the above two coefficient lists.
local f, b, i, zero;
    b := convert(evalm(linalg[multiply](linalg[matrix](bm),
        linalg[transpose](linalg[matrix](nullplus)))), listlist);
    zero := op(convert(linalg[matrix](1, nops(b[1]), 0), listlist));
    f := fminus;
    for i from nops(b) by -1 to 1 do
        if b[i] = zero then
            b := subsop(i = NULL, b); f := subsop(i = NULL, f)
        fi
    od;
    [f, b]
end;

connected := proc(G)
# Returns truw if graph is connected
local vert, vert1, vert2, vert3, x, y;
    if G = [] then RETURN(true) fi;
    vert := {op(map(op, G))};
    vert1 := {};
    vert2 := {vert[1]};
    while not (vert1 = vert2) do
        vert3 := {};
        for x in vert2 minus vert1 do for y in G do
                if member(x, y) then vert3 := vert3 union {op(y)} fi
            od
        od;
        vert1 := vert2;
        vert2 := vert2 union vert3
    od;
    evalb(vert2 = vert)
end;


nosplitsummand := proc(G)
#Returns true if there is no split summand of degree 2 or less
#for a graph created in possgraphs2 or if there is no split
#summand of degree 1 or no non-2-vertex split summand of 
#degree 2 in a graph created by possgraphs3.  See Lemma
#4.1 of [H-V].
local splitdeg, vert0, vert, vert1, vert2, vert3, x, y, H;
    vert0 :={op( vertices(G))};
    H := [];
    for x in G do if not member(1, x) then H := [op(H), x] fi od;
    vert := {op(vertices(H))};
    if nops(vert0) > nops(vert) + 1 then RETURN(false) fi; 
    splitdeg := 4;
    while (splitdeg > 2) and not(vert = {}) do
         vert1 := {};
         vert2 := {vert[1]};
         while not (vert1 = vert2) do
             vert3 := {};
             for x in vert2 minus vert1 do for y in H do
                     if member(x, y) then vert3 := vert3 union {op(y)} fi
                 od
             od;
             vert1 := vert2;
             vert2 := vert2 union vert3
         od; 
         if not(member(2,vert2)) then splitdeg := nops(vert2) fi;
         vert := vert minus vert2
    od;
    evalb(splitdeg > 2)
end;

calcul := proc(Gr, outfilename)
# Given list of graphs Gr, and the boundary map from deg-dimensional
# chains corresponding to maximal elements of degree deg to 
# (deg-1)-dimensional chains, this calculates the kernal of the
# boundary map.
local i, j, k, n, m, l, p, q, r, s, num, result, K, B, F, cubes, gr, gr2,
cols, BM, FM, nullsp, cubelist, kern, rows, divider, rowsp, colslist,
nonzero, rowsize, thismatrix, rowline, thisrow, kern2, cubelist2, tempk,
tempc;
    divider := `-----`;
    K := [];
    B := [];
    F := [];
    cubes := [];
    gr := [];
    gr2 := [];
    cols := 0;
    colslist := [];
    kern2 := [];
    cubelist2 := [];
    for num from nops(Gr) by -1 to 1 do
        appendto(outfilename);
        lprint(divider);
        lprint(`Examining element`, num, `of the list of possible graphs.`)
            ;
        lprint(`G is:`, Gr[num]);
        writeto(terminal);
        result := submain(Gr[num]);
        if result[1] = [] then
            subsop(num = NULL, Gr);
            appendto(outfilename);
            lprint(`Finished with graph`, num);
            writeto(terminal)
        elif result[2] = [] then
            appendto(outfilename);
            lprint(`Each element of kern listed below is in the kernel.`);
            lprint(`Elements of kern are indexed using the cubelist.`);
            lprint(`Kern is:`, result[1]);
            lprint(`Cubelist is:`, result[4]);
            lprint(`Finished with graph`, num);
            kern2 := [op(kern2), result[1]];
            cubelist2 := [op(cubelist2), result[4]];
            writeto(terminal);
            subsop(num = NULL, Gr)
        else
            nullsp := linalg[nullspace](linalg[matrix](result[2]));
            if nops(nullsp) = 0 then
                K := [op(K), result[1]];
                B := [op(B), result[2]];
                F := [op(F), result[3]];
                cubes := [op(cubes), result[4]];
                cols := cols + nops(result[2][1]);
                colslist := [op(colslist), nops(result[2][1])];
                appendto(outfilename);
                lprint(`Saving the following for further computations:`);
                lprint(`K is:`, result[1]);
                lprint(`B is:`, result[2]);
                lprint(`F is:`, result[3]);
                lprint(`Cubes are:`, result[4]);
                writeto(terminal)
            else
                nullsp := linalg[matrix](
                    [op(map(v -> convert(v, list), nullsp))]);
                kern := convert(
                    linalg[multiply](nullsp, linalg[matrix](result[1])),
                    listlist);
                tempk := kern;
                tempc := result[4];
                for q from nops(tempk[1]) by -1 to 1 do
                    nonzero := false;
                    for r to nops(tempk) do
                        if not (tempk[r][q] = 0) then nonzero := true fi
                    od;
                    if nonzero = false then
                        for r to nops(tempk) do tempk :=
                            subsop(r = subsop(q = NULL, tempk[r]), tempk)
                        od;
                        tempc := subsop(q = NULL, tempc)
                    fi
                od;
                kern2 := [op(kern2), tempk];
                cubelist2 := [op(cubelist2), tempc];
                appendto(outfilename);
                lprint(
                    `Each element of kern listed below is in the kernel.`)
                    ;
                lprint(`Elements of kern are indexed using the cubelist.`)
                    ;
                lprint(`Kern is:`, tempk);
                lprint(`Cubelist is:`, tempc);
                writeto(terminal);
                rowsp := linalg[rowspace](linalg[matrix](result[2]));
                if nops(rowsp) = 0 then subsop(num = NULL, Gr)
                else
                    cols := cols + nops(rowsp);
                    colslist := [op(colslist), nops(rowsp)];
                    rowsp := linalg[matrix](
                        [op(map(v -> convert(v, list), rowsp))]);
                    K := [op(K), convert(
                        linalg[multiply](rowsp, linalg[matrix](result[1]))
                        , listlist)];
                    B := [op(B), convert(linalg[multiply](
                        linalg[matrix](result[2]),
                        linalg[transpose](rowsp)), listlist)];
                    F := [op(F), result[3]];
                    cubes := [op(cubes), result[4]];
                    n := nops(K);
                    appendto(outfilename);
                    lprint(
                        `Saving the following for further computations:`);
                    lprint(`K is:`, K[n]);
                    lprint(`B is:`, B[n]);
                    lprint(`F is:`, F[n]);
                    lprint(`Cubes are:`, cubes[n]);
                    writeto(terminal)
                fi
            fi
        fi
    od;
    if 0 < cols then
        appendto(outfilename);
        lprint(divider);
        lprint(`Beginning to combine previously saved results.`);
        writeto(terminal);
        gc();
        BM := [];
        rowsize := 0;
        for l to nops(K) do
            appendto(outfilename);
            lprint(`Echelon form calculations.  Stage:`, l);
            writeto(terminal);
            thismatrix := [];
            for i to nops(B[l]) do
                rowline := F[l][i];
                thisrow := B[l][i];
                for j from l + 1 to nops(K) do
                    if member(rowline, F[j], 'm') then
                        F := subsop(j = subsop(m = NULL, F[j]), F);
                        thisrow := [op(thisrow), op(B[j][m])];
                        B := subsop(j = subsop(m = NULL, B[j]), B)
                    else thisrow :=
                        [op(thisrow), op(zerolist(colslist[j]))]
                    fi
                od;
                thismatrix := [op(thismatrix), thisrow]
            od;
            if not (thismatrix = []) then
                thismatrix := linalg[rowspace](linalg[matrix](thismatrix))
                    ;
                thismatrix := [op(map((M, rowsize) ->
                    [op(zerolist(rowsize)), op(convert(M, list))],
                    thismatrix, rowsize))]
            fi;
            rowsize := rowsize + colslist[l];
            BM := [op(BM), op(thismatrix)]
        od;
        appendto(outfilename);
        lprint(`Done with Echelon form calculations.  Calculating nullsp.`)
            ;
        writeto(terminal);
        nullsp := linalg[nullspace](linalg[matrix](BM));
        nullsp := map(M -> convert(M, list), nullsp);
        appendto(outfilename);
        lprint(`Nullsp calculated.`);
        lprint(`Nullsp is:`, nullsp);
        lprint(`Using nullsp to try to find elements in the kernel.`);
        writeto(terminal);
        cubelist := sort([op({op(map(op, cubes))})], lex);
        rows := nops(nullsp);
        if rows = 0 then
            appendto(outfilename);
            lprint(`No more elements found in the kernel.`);
            writeto(terminal);
            RETURN([kern2, cubelist2])
        fi;
        kern := linalg[matrix](rows, nops(cubelist), 0);
        cols := 0;
        for l to nops(K) do
            for i to nops(cubes[l]) do
                member(cubes[l][i], cubelist, 'm');
                for j to nops(K[l]) do for n to nops(nullsp) do kern[n, m]
                         := kern[n, m] + nullsp[n][cols + j]*K[l][j][i]
                    od
                od
            od;
            cols := cols + colslist[l]
        od;
        kern := convert(kern, listlist);
        for q from nops(cubelist) by -1 to 1 do
            nonzero := false;
            for r to nops(kern) do
                if not (kern[r][q] = 0) then nonzero := true fi
            od;
            if nonzero = false then
                for r to nops(kern) do
                    kern := subsop(r = subsop(q = NULL, kern[r]), kern)
                od;
                cubelist := subsop(q = NULL, cubelist)
            fi
        od;
        kern2 := [op(kern2), kern];
        cubelist2 := [op(cubelist2), cubelist];
        appendto(outfilename);
        lprint(divider);
        lprint(`Each element of kern listed below is in the kernel.`);
        lprint(`Elements of kern are indexed using the cubelist.`);
        lprint(`Kern is:`, kern);
        lprint(`Cubelist is:`, cubelist);
        writeto(terminal)
    fi;
    [kern2, cubelist2]
end;

cleanitup := proc(result, cubesused)
# Writes out the kernel kern for main by properly reindexing
# the results of calcul
local i, j, k, m, kern, elekern;
    kern := [];
    for i to nops(result[1]) do for j to nops(result[1][i]) do
            elekern := [];
            for k to nops(cubesused) do
                if member(cubesused[k], result[2][i], 'm') then
                    elekern := [op(elekern), result[1][i][j][m]]
                else elekern := [op(elekern), 0]
                fi
            od;
            kern := [op(kern), elekern]
        od
    od;
    kern
end;

                    supsort := L -> sort(insort(L), lexico);

                comb := (L, k) -> combinat[choose](edges(L), k);

core := proc(L)
# Keep removing free edges inductively until left with core graph
local G, H;
    H := L; G := trim(L); while not (G = H) do H := G; G := trim(G) od; H
end;

                            e := G -> nops({op(G)});

                       v := G -> nops({op(map(op, G))});

findfplusandbp := proc(plusfaces, plusfolds, bplus)
# Remove plusfolds from plusface list and clean it up
local fplus, bp, i, zero;
  if plusfaces = [] then
     [[],[]]
  else
    fplus := plusfaces;
    bp := convert(bplus, listlist);
    zero := op(convert(linalg[matrix](1, nops(bp[1]), 0), listlist));
    for i from nops(bp) by -1 to 1 do
        if bp[i] = zero then
            fplus := subsop(i = NULL, fplus); bp := subsop(i = NULL, bp)
        fi
    od;
    for i from nops(bp) by -1 to 1 do
        if member(fplus[i], plusfolds) then
            fplus := subsop(i = NULL, fplus); bp := subsop(i = NULL, bp)
        fi
    od;
    [fplus, bp]
  fi
end;

addones := proc(G)
# Used in possgraphs2 and possgraphs3 which start out by
# omitting the 1's in graph lists
local x, H;
    H := [];
    for x in G do
        if not type(x, list) then H := [op(H), [1, x]]
        else H := [op(H), x]
        fi
    od;
    H
end;

findbplus := proc(vert, P, cubes, plusfaces)
# Create matrix listing plusface boundaries of cubes
local bplus, X, i, m, n;
    bplus := linalg[matrix](nops(plusfaces), nops(cubes), 0);
    for i to vert - 1 do for n to nops(cubes) do
            X := bdryplus(cubes[n], i, P);
            member(X[1 .. 2], plusfaces, 'm');
            bplus[m, n] := bplus[m, n] + X[3]
        od
    od;
    bplus
end;

                   edges := G -> sort([op({op(G)})], lexico);

findnullplus := proc(bp)
# Find nullspace of listlist bp, value returned is listlist
    sort([
    op(map(v -> convert(v, list), linalg[nullspace](linalg[matrix](bp))))]
    , lex)
end;

calcul2 := proc(Gr, cubesused, outfilename)
# Uses the degree deg+1 graphs in Gr to find higher dimensional
# chains of cubes whose boundaries are entirely contained in
# cubesused, and gives a basis for the image of such chains
local i, j, k, n, m, l, p, q, r, s, t, num, result, K, B, F, cubes, gr,
gr2, cols, BM, FM, nullsp, cubelist, kern, rows, divider, rowsp, colslist,
bcut, elebcut, imlist, rowsize, thismatrix, rowline, thisrow, colsp,
kerfin, imfin, nonzero, imaug, zero;
    divider := `-----`;
    imfin := [];
    kerfin := [];
    cubelist := [];
    K := [];
    B := [];
    F := [];
    cubes := [];
    gr := [];
    gr2 := [];
    cols := 0;
    colslist := [];
    bcut := [];
    for num from nops(Gr) by -1 to 1 do
        appendto(outfilename);
        lprint(divider);
        lprint(`Examining element`, num,
            `of the list of higher degree graphs.`);
        lprint(`G is:`, Gr[num]);
        writeto(terminal);
        result := submain(Gr[num]);
        if result[1] = [] then
            subsop(num = NULL, Gr);
            appendto(outfilename);
            lprint(`Finished with graph`, num);
            writeto(terminal)
        elif result[2] = [] then
            appendto(outfilename);
            lprint(`Finished with graph`, num);
            writeto(terminal);
            subsop(num = NULL, Gr)
        else
            nullsp := linalg[nullspace](linalg[matrix](result[2]));
            if nops(nullsp) = 0 then
                for q to nops(result[1]) do
                    elebcut := [];
                    for r to nops(cubesused) do
                        if member(cubesused[r], result[3], 'm') then
                            elebcut := [op(elebcut), result[2][m][q]];
                            if q = nops(result[1]) then result := subsop(
                                2 = subsop(m = NULL, result[2]),
                                3 = subsop(m = NULL, result[3]), result)
                            fi
                        else elebcut := [op(elebcut), 0]
                        fi
                    od;
                    bcut := [op(bcut), elebcut]
                od;
                if nops(result[2]) = 0 then result := subsop(
                    2 = [zerolist(nops(result[1]))], 3 = [spacer], result)
                fi;
                K := [op(K), result[1]];
                B := [op(B), result[2]];
                F := [op(F), result[3]];
                cubes := [op(cubes), result[4]];
                cols := cols + nops(result[2][1]);
                colslist := [op(colslist), nops(result[2][1])];
                appendto(outfilename);
                lprint(`Saving the following for further computations:`);
                lprint(`K is:`, result[1]);
                lprint(`B is:`, result[2]);
                lprint(`F is:`, result[3]);
                lprint(`Cubes are:`, result[4]);
                writeto(terminal)
            else
                nullsp := linalg[matrix](
                    [op(map(v -> convert(v, list), nullsp))]);
                kern := convert(
                    linalg[multiply](nullsp, linalg[matrix](result[1])),
                    listlist);
                rowsp := linalg[rowspace](linalg[matrix](result[2]));
                if nops(rowsp) = 0 then subsop(num = NULL, Gr)
                else
                    cols := cols + nops(rowsp);
                    colslist := [op(colslist), nops(rowsp)];
                    rowsp := linalg[matrix](
                        [op(map(v -> convert(v, list), rowsp))]);
                    K := [op(K), convert(
                        linalg[multiply](rowsp, linalg[matrix](result[1]))
                        , listlist)];
                    B := [op(B), convert(linalg[multiply](
                        linalg[matrix](result[2]),
                        linalg[transpose](rowsp)), listlist)];
                    F := [op(F), result[3]];
                    cubes := [op(cubes), result[4]];
                    n := nops(K);
                    for q to nops(K[n]) do
                        elebcut := [];
                        for r to nops(cubesused) do
                            if member(cubesused[r], F[n], 'm') then
                                elebcut := [op(elebcut), B[n][m][q]];
                                if q = nops(K[n]) then
                                    B := subsop(
                                        n = subsop(m = NULL, B[n]), B);
                                    F := subsop(
                                        n = subsop(m = NULL, F[n]), F)
                                fi
                            else elebcut := [op(elebcut), 0]
                            fi
                        od;
                        bcut := [op(bcut), elebcut]
                    od;
                    if nops(B[n]) = 0 then
                        B := subsop(n = [zerolist(nops(K[n]))], B);
                        F := subsop(n = [spacer], F)
                    fi;
                    appendto(outfilename);
                    lprint(
                        `Saving the following for further computations:`);
                    lprint(`K is:`, K[n]);
                    lprint(`B is:`, B[n]);
                    lprint(`F is:`, F[n]);
                    lprint(`Cubes are:`, cubes[n]);
                    writeto(terminal)
                fi
            fi
        fi
    od;
    if 0 < cols then
        appendto(outfilename);
        lprint(divider);
        lprint(`Beginning to combine previously saved results.`);
        lprint(`Bcut is:`, bcut);
        writeto(terminal);
        gc();
        BM := [];
        rowsize := 0;
        for l to nops(K) do
            appendto(outfilename);
            lprint(`Echelon form calculations.  Stage:`, l);
            writeto(terminal);
            thismatrix := [];
            for i to nops(B[l]) do
                rowline := F[l][i];
                thisrow := B[l][i];
                for j from l + 1 to nops(K) do
                    if member(rowline, F[j], 'm') then
                        F := subsop(j = subsop(m = NULL, F[j]), F);
                        thisrow := [op(thisrow), op(B[j][m])];
                        B := subsop(j = subsop(m = NULL, B[j]), B)
                    else thisrow :=
                        [op(thisrow), op(zerolist(colslist[j]))]
                    fi
                od;
                thismatrix := [op(thismatrix), thisrow]
            od;
            if not (thismatrix = []) then
                thismatrix := linalg[rowspace](linalg[matrix](thismatrix))
                    ;
                thismatrix := [op(map((M, rowsize) ->
                    [op(zerolist(rowsize)), op(convert(M, list))],
                    thismatrix, rowsize))]
            fi;
            rowsize := rowsize + colslist[l];
            BM := [op(BM), op(thismatrix)]
        od;
        appendto(outfilename);
        lprint(`Done with Echelon form calculations.  Calculating nullsp.`)
            ;
        writeto(terminal);
        nullsp := linalg[nullspace](linalg[matrix](BM));
        nullsp := map(M -> convert(M, list), nullsp);
        appendto(outfilename);
        lprint(`Nullsp calculated.`);
        lprint(`Nullsp is:`, nullsp);
        lprint(`Using nullsp to try to find elements in the image.`);
        writeto(terminal);
        if nops(nullsp) = 0 then
            appendto(outfilename);
            lprint(`No elements found in the image.`);
            writeto(terminal);
            RETURN([[], [], []])
        fi;
        cubelist := sort([op({op(map(op, cubes))})], lex);
        rows := nops(nullsp);
        kern := linalg[matrix](rows, nops(cubelist), 0);
        cols := 0;
        for l to nops(K) do
            for i to nops(cubes[l]) do
                member(cubes[l][i], cubelist, 'm');
                for j to nops(K[l]) do for n to nops(nullsp) do kern[n, m]
                         := kern[n, m] + nullsp[n][cols + j]*K[l][j][i]
                    od
                od
            od;
            cols := cols + colslist[l]
        od;
        nullsp := linalg[matrix]([op(nullsp)]);
        bcut := linalg[matrix](bcut);
        imlist := linalg[multiply](nullsp, bcut);
        imaug :=
            linalg[rref](linalg[augment](imlist, kern), nops(cubesused));
        imfin := linalg[delcols](imaug,
            nops(cubesused) + 1 .. nops(cubesused) + nops(cubelist));
        kerfin := linalg[delcols](imaug, 1 .. nops(cubesused));
        imfin := convert(imfin, listlist);
        kerfin := convert(kerfin, listlist);
        zero := zerolist(nops(cubesused));
        t := nops(imfin);
        for q from t by -1 to 1 do
            if imfin[q] = zero then
                imfin := subsop(q = NULL, imfin);
                kerfin := subsop(q = NULL, kerfin)
            fi
        od;
        for q from nops(cubelist) by -1 to 1 do
            nonzero := false;
            for r to nops(kerfin) do
                if not (kerfin[r][q] = 0) then nonzero := true fi
            od;
            if nonzero = false then
                for r to nops(kerfin) do kerfin :=
                    subsop(r = subsop(q = NULL, kerfin[r]), kerfin)
                od;
                cubelist := subsop(q = NULL, cubelist)
            fi
        od;
        appendto(outfilename);
        lprint(divider);
        lprint(`Each element of imfin listed below is in the image.`);
        lprint(`Elements of imfin are indexed using cubesused.`);
        lprint(`Each element of kerfin listed below is a chain of`);
        lprint(`degree deg+1 that maps onto the corresponding`);
        lprint(`element of imfin.`);
        lprint(`Elements of kerfin are indexed using cubelist.`);
        lprint(`Imfin is:`, imfin);
        lprint(`Cubesused is:`, cubesused);
        lprint(`Kerfin is:`, kerfin);
        lprint(`Cubelist is:`, cubelist);
        writeto(terminal)
    fi;
    [imfin, kerfin, cubelist]
end;


nperm := proc(X)
# Returns nice permutations which when applied to X give the
# normal form of X (after possibly being supsorted)
local P, Y;
    P := perm(X);
    Y := nf(X);
    select((M, X, Y) -> supsort(cd(M, X)) = Y, P, X, Y)
end;

          findnqgraphs := qr -> sort([op({op(map(nf, qr[1]))})], lex);
# Looks at graphs involved in qr, takes their normal forms and puts
# them into a sorted list with no duplications

main := proc(deg, outfilename, finalfilename)
# The main procedure which calculates the stable rational homology of the
# holomorph F_n \rtimes Aut(F_n) in dimension deg; intermediate results
# are written to outfilename, and the final results are written to
# finalfilename
local result, result2, Gr, Gr2, divider, kern, cubesused;
    if type(deg, integer) and 1 < deg and deg < 7 then
        divider := `-----`;
        writeto(outfilename);
        lprint(`deg is:`, deg);
        lprint(divider);
        lprint(`Beginning calculation of possible graphs.`);
        writeto(terminal);
        Gr := possgraphs(deg, outfilename);
        appendto(outfilename);
        lprint(divider);
        lprint(`Possible graphs are:`, Gr);
        lprint(divider);
        lprint(`Beginning calculation of kernel.`);
        writeto(terminal);
        result := calcul(Gr, outfilename);
        if not (result[2] = []) then
            cubesused := sort([op({op(map(op, result[2]))})], lex);
            kern := cleanitup(result, cubesused);
            appendto(outfilename);
            lprint(divider);
            lprint(`Starting with preliminaries of calculating the`);
            lprint(`image of the higher degree graphs.`);
            lprint(`Beginning calculation of higher degree graphs.`);
            writeto(terminal);
            Gr2 := possgraphs(deg + 1, outfilename);
            appendto(outfilename);
            lprint(divider);
            lprint(`Higher degree graphs are:`, Gr2);
            lprint(divider);
            lprint(`Beginning calculation of image.`);
            writeto(terminal);
            result2 := calcul2(Gr2, cubesused, outfilename)
        fi;
        appendto(outfilename);
        lprint(divider);
        lprint(`Calculations completed.`);
        writeto(terminal);
        if finalfilename = outfilename then
            appendto(outfilename); lprint(divider); lprint(divider)
        else writeto(finalfilename)
        fi;
        lprint(`Final statement of results:`);
        lprint(`Deg is:`, deg);
        if result[2] = [] then
            lprint(`The homology group H_deg of the quotient of auter`);
            lprint(`space is zero.`)
        else
            lprint(
                `Elements of kern listed below span the top dimensional`);
            lprint(`homology group of the quotient complex Q_n,deg for \
                n large.`);
            lprint(`Elements of kern are indexed by cubesused.`);
            lprint(`Kern is:`, kern);
            lprint(`Cubesused is:`, cubesused);
            lprint(`We now calculate the image imfin of deg+1 chains in`);
            lprint(
                `Q_n,deg+1 that map to kern.  Elements of imfin are also`)
                ;
            lprint(`indexed by cubesused.`);
            lprint(
                `For each element of imfin, we give a corresponding deg+1`)
                ;
            lprint(
                `chain that maps onto it.  The deg+1 chains are listed in`)
                ;
            lprint(
                `kerfin.  Elements of kerfin are indexed by "cubelist".`);
            lprint(`Imfin is:`, result2[1]);
            lprint(`Kerfin is:`, result2[2]);
            lprint(`Cubelist is:`, result2[3])
        fi;
        writeto(terminal)
    else
        lprint(`ERROR.  First argument must be an integer greater than 1 `)
            ;
        lprint(`and less than 7.`);
        lprint(
            `Second and third arguments must be a filenames for output,`);
        lprint(`hence either strings enclosed in backquotes or terminal.`)
    fi
end;


numboc := proc(L, M)
# The number of times M occurs in L
local i, j;
    i := 0; for j to nops(L) do if L[j] = M then i := i + 1 fi od; i
end;

                         forestq := L -> core(L) = [];

findbminus := proc(vert, cubes, qr, minusfaces)
# Returns a minusfaces x cubes dimensional matrix
# listing the minusface boundaries of cubes
local bminus, X, i, m, n;
    bminus := linalg[matrix](nops(minusfaces), nops(cubes), 0);
    for i to vert - 1 do for n to nops(cubes) do
            X := byminus(cubes[n], i, qr);
            member(X[1 .. 2], minusfaces, 'm');
            bminus[m, n] := bminus[m, n] + X[3]
        od
    od;
    bminus
end;

reflq := proc(Y, P)
# Is there a reflection in P inducing an odd permutation of the edges of Y?
local i, Z, compare1, compare2;
    compare1 := supsort(Y);
    compare2 := -signature(insort(Y));
    for i to nops(P) do
        Z := map((M, i, P) -> map((k, i, P) -> P[i][k], M, i, P), Y, i, P)
            ;
        if supsort(Z) = compare1 and signature(insort(Z)) = compare2 then
            RETURN(true)
        fi
    od;
    false
end;

findcubes := proc(G, trees)
# See Lemma 4.3 and Lemma 4.4 of [H-V]
# With graph G and tree list, forms list of pairs consisting
# of the graph and the tree.  From these cubes, we select
# out those which satisfy nobasetriangleq and notriangleq
local cubes, i;
    cubes := [];
    for i to nops(trees) do cubes := [op(cubes), [G, trees[i]]] od;
    cubes := select(nobasetriangleq, select(notriangleq, cubes));
    cubes
end;

notfreeq :=

    (L, M) -> 1 < numboc(map(op, L), M[1]) and 1 < numboc(map(op, L), M[2]);
                         


aut := proc(X)
local P;
    P := perm(X); select((M, X) -> supsort(cd(M, X)) = supsort(X), P, X)
end;

perm := proc(G)
local X, satisfied;
    X := supsort(G);
    satisfied := L -> L[1] = 1;
    ok(sigma(max(op(vertices(G)))), {2, op(map(op, select(satisfied, X)))})
end;

mpnf := proc(X, P)
local i, Z, W, a, graphset;
    Z := supsort(X);
    graphset := {};
    for i to nops(P) do
        a := (M, i, P) -> map((k, i, P) -> P[i][k], M, i, P);
        W := supsort(map(a, X, i, P));
        graphset := graphset union {W};
        if lex(W, Z) then Z := W fi
    od;
    [Z, graphset]
end;

lex := proc(l1, l2)
    if not (type(l1, list) or type(l2, list)) then RETURN(evalb(l1 <= l2))
    elif not type(l1, list) and type(l2, list) then RETURN(true)
    elif type(l1, list) and not type(l2, list) then RETURN(false)
    elif l1 = [] then RETURN(true)
    elif l2 = [] then RETURN(false)
    elif lex(l1[1], l2[1]) then
        if not (l1[1] = l2[1]) then RETURN(true)
        elif nops(l1) = 1 then RETURN(true)
        elif nops(l2) = 1 then RETURN(false)
        else RETURN(lex(l1[2 .. nops(l1)], l2[2 .. nops(l2)]))
        fi
    else RETURN(false)
    fi
end;

qrcreate := proc(G)
local Q, R, i, edgelist;
    Q := [];
    R := [];
    edgelist := edges(G);
    for i to e(G) do
        Q := [op(Q), collapse(G, edgelist[i])]; R := [op(R), nperm(Q[i])]
    od;
    [Q, R]
end;

        cd := (X, G) -> map((M, X) -> map((k, X) -> X[k], M, X), G, X);
#cd applies perm X to graph G

ad := proc(L1, L2)
local i, M;
    M := []; for i to nops(L1) do M := [op(M), L1[i] + L2[i]] od; M
end;

            findfolds := (P, forests) -> select(reflq, forests, P);

mmu := proc(L1, L2)
local i, j, k, M, N, x;
    M := [];
    for i to nops(L1) do
        N := [];
        for j to nops(L2[1]) do
            x := 0;
            for k to nops(L2) do x := x + L1[i][k]*L2[k][j] od;
            N := [op(N), x]
        od;
        M := [op(M), N]
    od;
    M
end;

sigma := proc(n)
local one;
    one := L -> [1, 2, op(L)]; map(one, combinat[permute]([`$`(3 .. n)]))
end;

pnf := proc(X, P)
local i, Z, W, a;
    Z := supsort(X);
    for i to nops(P) do
        a := (M, i, P) -> map((k, i, P) -> P[i][k], M, i, P);
        W := supsort(map(a, X, i, P));
        if lex(W, Z) then Z := W fi
    od;
    Z
end;

                      nf := X -> pnf(X, perm(supsort(X)));


byminus := proc(L, k, QR)
local X, Y, p, P, R, Q;
    X := L[1];
    Y := L[2];
    Q := QR[1];
    R := QR[2];
    member(Y[k], edges(X), 'p');
    P := R[p];
    pnormalform([Q[p], collapse(Y, Y[k]), (-1)^k], P)
end;

ok := proc(P, S)
local v, X, satisfied;
    X := P;
    satisfied := (L, S, v) -> L[S[v]] <= nops(S);
    for v from 2 to nops(S) do X := select(satisfied, X, S, v) od;
    X
end;

nobasetriangleq := proc(L)
# Check for a triangle in the graph L[1] at the basepoint
# where the subforest L[2] contains exactly one edge of this 
# triangle, and where the triangle has a double edge adjacent
# to the basepoint, and where the edge that the forest contains
# is not the single edge connected to the basepoint
# See Lemma 4.4 of [H-V]
local n, i, j, G, F, selectfn;
    G := L[1];
    F := L[2];
    n := v(G);
    selectfn := (M, i) -> evalb(M = [1, i]);
    for i from 2 to n do for j from i + 1 to n do
            if nops(select(selectfn, G, i)) = 2 and
            nops({[i, j], [1, i], [1, j]} intersect {op(G)}) = 3 and (
            {[i, j], [1, i], [1, j]} intersect {op(F)} = {[1, i]} or
            {[i, j], [1, i], [1, j]} intersect {op(F)} = {[i, j]}) then
                RETURN(false)
            fi
        od
    od;
    true
end;

mu := proc(a, L)
# Multiply vector L by scalar m
local i, M;
    M := []; for i to nops(L) do M := [op(M), a*L[i]] od; M
end;


findplusfolds := proc(G, folds)
# Find plusfaces that are folds and thus can be ignored
local plusfolds, i;
    plusfolds := [];
    for i to nops(folds) do plusfolds := [op(plusfolds), [G, folds[i]]] od
    ;
    plusfolds
end;

findtrees := proc(G, vert, P)
# Find isomorphism classes of trees in G
local L;
    L := select(noreflq,{op(map(pnf, 
         select(forestq, comb(G, vert - 1)), P))},P);
    sort([op(L)], lex)
end;

              trim := L -> select((M, L) -> notfreeq(L, M), L, L);
# Trim off free edges from graph

findfminusandbm := proc(minusfaces, minusfolds, bminus)
# Nix minusfolds from list of minusfaces
local fminus, bm, i, zero;
    fminus := minusfaces;
    bm := convert(bminus, listlist);
    zero := op(convert(linalg[matrix](1, nops(bm[1]), 0), listlist));
    for i from nops(bm) by -1 to 1 do
        if bm[i] = zero then
            fminus := subsop(i = NULL, fminus); bm := subsop(i = NULL, bm)
        fi
    od;
    for i from nops(bm) by -1 to 1 do
        if member(fminus[i], minusfolds) then
            fminus := subsop(i = NULL, fminus); bm := subsop(i = NULL, bm)
        fi
    od;
    [fminus, bm]
end;

findforests := proc(G, vert, P, trees)
# Find forests, which give plusfaces
local K;
    K := {op(map(op, map(combinat[choose], trees, vert - 2)))};
    K := sort([op(map(pnf, K, P))], lex);
    K
end;

possgraphs := proc(deg, outfilename)
# maximal degree deg graphs
    select(nosplitsummand,  
          possgraphs2(deg,outfilename) union possgraphs3(deg-1,outfilename));
end;

possgraphs2 := proc(deg, outfilename)
# maximal degree deg graphs with vertex 1 = vertex 2
local L, K, i, j, n, G, G1, H, inserts, slots, x, Ldone;
     L :=    {[3, 4, 5, [3, 5], [4, 5]], 
[3, 4, 4, 5, 5, [3, 5]], 
[3, 3, 4, 4, 5, 5, 5], 
[3, 3, 4, 5, 5, [4, 5]], 
[3, 4, 5, 5, 5, [3, 4]], 
[4, 5, 5, [3, 4], [3, 5]], 
[3, 5, 5, [3, 4], [4, 5]], 
[5, [3, 4], [3, 5], [4, 5]]};
    for i from 6 to deg + 2 do
        appendto(outfilename);
        lprint(`Adding vertex`, i, `to possible graphs.`);
        writeto(terminal);
        K := {};
        for G in L do
            slots := [];
            Ldone := [];
            for n to nops(G) do
                if not type(G[n], list) and not member(G[n], Ldone) then 
                    slots := [op(slots), n]; 
                    Ldone := [op(Ldone), G[n]]
                fi
            od;
            H := [[op(G), i, i, i]];
            inserts := combinat[choose](slots, 1);
            if inserts = [[]] then inserts := [] fi;
            for x in inserts do
                H := [op(H), [op(subsop(x[1] = [G[x[1]], i], G)), i, i]]
            od;
            inserts := combinat[choose](slots, 2);
            if inserts = [[]] then inserts := [] fi;
            for x in inserts do
              G1 := [op(subsop(x[1] = [G[x[1]], i], x[2] = [G[x[2]], i], G)),
                    i];
                H := [op(H), G1]
            od;
            inserts := combinat[choose](slots, 3);
            if inserts = [[]] then inserts := [] fi;
            for x in inserts do
                G1 := subsop(x[1] = [G[x[1]], i], x[2] = [G[x[2]], i],
                    x[3] = [G[x[3]], i], G);
                H := [op(H), G1]
            od;
            H := {op(map(x -> sort(x, lex), H))};
            K := K union H
        od;
        L := K
    od;
    L := map(add34, L);
    appendto(outfilename);
    lprint(`Adding 1's.`);
    writeto(terminal);
    L := map(addones, L);
    L := {op(map(x -> sort(x, lex), L))};
    appendto(outfilename);
    lprint(`Removing graphs which have a separating edge.`);
    writeto(terminal);
    L := select(nosepedge, L);
    appendto(outfilename);
    lprint(`Putting graphs in normal form.`);
    writeto(terminal);
    L := map(G -> op(zapgraph2(G)), L);
    L
end;

possgraphs3 := proc(deg, outfilename)
# maximal degree deg + 1 graphs with vertex 1 not equal
# to vertex 2 and vertex 2 having valence 2
local L, K, i, j, n, G, G1, H, inserts, slots, x, Ldone;
    L :=    {[4, [2, 3], [2, 4], [3, 4]], 
[2, 3, 4, 4, 4, [2, 3]], 
[3, 4, 4, [2, 3], [2, 4]], 
[2, 2, 3, 3, 4, 4, 4], 
[2, 4, 4, [2, 3], [3, 4]], 
[2, 3, 3, 4, 4, [2, 4]], 
[2, 2, 3, 4, 4, [3, 4]], 
[2, 3, 4, [2, 4], [3, 4]]};
    for i from 5 to deg + 2 do
        appendto(outfilename);
        lprint(`Adding vertex`, i, `to possible graphs.`);
        writeto(terminal);
        K := {};
        for G in L do
            slots := [];
            Ldone := [];
            for n to nops(G) do
                if not type(G[n], list) and not member(G[n], Ldone) then 
                    slots := [op(slots), n]; 
                    Ldone := [op(Ldone), G[n]]
                fi
            od;  
            H := [[op(G), i, i, i]];
            inserts := combinat[choose](slots, 1);
            if inserts = [[]] then inserts := [] fi;
            for x in inserts do
                H := [op(H), [op(subsop(x[1] = [G[x[1]], i], G)), i, i]]
            od;
            inserts := combinat[choose](slots, 2);
            if inserts = [[]] then inserts := [] fi;
            for x in inserts do
                G1 := [op(
                    subsop(x[1] = [G[x[1]], i], x[2] = [G[x[2]], i], G)),
                    i];
                H := [op(H), G1]
            od;
            inserts := combinat[choose](slots, 3);
            if inserts = [[]] then inserts := [] fi;
            for x in inserts do
                G1 := subsop(x[1] = [G[x[1]], i], x[2] = [G[x[2]], i],
                    x[3] = [G[x[3]], i], G);
                H := [op(H), G1]
            od;
            H := {op(map(x -> sort(x, lex), H))};
            K := K union H
        od;
        L := K
    od;
    L := map(add3, L);
    appendto(outfilename);
    lprint(`Adding 1's.`);
    writeto(terminal);
    L := map(addones, L);
    L := {op(map(x -> sort(x, lex), L))};
     appendto(outfilename);
    lprint(`Removing graphs where the basepoint has valence < 2.`);
    writeto(terminal);
    L := select(G -> 1 < numboc(map(op, G), 1), L);
    appendto(outfilename);
    lprint(`Removing graphs which have a separating edge.`);
    writeto(terminal);
    L := select(nosepedge, L);
    appendto(outfilename);
    lprint(`Putting graphs in normal form.`);
    writeto(terminal);
    L := map(G -> op(zapgraph2(G)), L);
    L
end;


                    vertices := G -> [op({op(map(op, G))})];

add34 := proc(G)
    [3, 4, op(G)]
end;

add3 := proc(G)
    [3, op(G)]
end;

submain := proc(H)
# Called by both calcul and calcul2
# Looks at a particular graph H and finds cubes, plusfaces, minusfaces, etc,
# for it, as well as finding a basis for chains whose boundary evaluates
# to zero on the plusfaces.
local vert, P, G, trees, cubes, forests, plusfaces, folds, plusfolds,
bplus, fplusandbp, fplus, bp, nullplus, qr, nqgraphs, minusfacesandfolds,
minusfaces, minusfolds, bminus, fminusandbm, fminus, bm, fandb, f, b;
    vert := v(H);
#    G := nf(H); #Possgraphs already puts in normal form
    G := H;
    P := aut(G);
    trees := findtrees(G, vert, P);
    cubes := findcubes(G, trees);
    forests := findforests(G, vert, P, trees);
    plusfaces := findplusfaces(G, forests);
    folds := findfolds(P, forests);
    plusfolds := findplusfolds(G, folds);
    bplus := evalm(findbplus(vert, P, cubes, plusfaces));
    fplusandbp := findfplusandbp(plusfaces, plusfolds, bplus);
    fplus := fplusandbp[1];
    bp := fplusandbp[2];
    nullplus := findnullplus(bp);
    if nops(nullplus) = 0 then RETURN([[], [], [], cubes])
    else
        qr := qrcreate(G);
        nqgraphs := findnqgraphs(qr);
        minusfacesandfolds := findminusfacesandfolds(vert, nqgraphs);
        minusfaces := minusfacesandfolds[1];
        minusfolds := minusfacesandfolds[2];
        bminus := evalm(findbminus(vert, cubes, qr, minusfaces));
        fminusandbm := findfminusandbm(minusfaces, minusfolds, bminus);
        fminus := fminusandbm[1];
        bm := fminusandbm[2];
        fandb := findfandb(nullplus, fminus, bm);
        f := fandb[1];
        b := fandb[2];
        RETURN([nullplus, b, f, cubes])
    fi
end;

findminusfacesandfolds := proc(vert, nqgraphs)
# Finds minusfaces (faces ontained by collapsing an edge from the larger
# cube that we have to worry about) and minusfolds (faces obtained by
# collapsing an edge from the larger cube that we don't have to worry
# about since they are "folded over" and glued to themselves.)
local i, j, X, P, L, trees, folds, minusfaces, minusfolds;
    minusfaces := [];
    minusfolds := [];
    for i to nops(nqgraphs) do
        X := nqgraphs[i];
        P := aut(X);
        L := sort([op({op(map(pnf, comb(X, vert - 2), P))})], lex);
        trees := select(forestq, L);
        folds := select(reflq, trees, P);
        for j to nops(trees) do
            minusfaces := [op(minusfaces), [X, trees[j]]]
        od;
        for j to nops(folds) do
            minusfolds := [op(minusfolds), [X, folds[j]]]
        od
    od;
    [minusfaces, minusfolds]
end;

nosepedge := proc(G)
# Does the graph have a separating edge?
local i;
    for i to nops(G) do
        if not connected(subsop(i = NULL, G)) then RETURN(false) fi
    od;
    true
end;

removebase := proc(G)
# Remove basepoint from graph and associated edges
local x, H, vert, vert1;
    vert := vertices(G);
    H := [];
    for x in G do if not member(1, x) then H := [op(H), x] fi od;
    vert1 := [1, op(vertices(H))];
    for x in {op(vert)} minus {op(vert1)} do H := [op(H), [x, x]] od;
    H
end;

                          insort := L -> map(sort, L);

                     noreflq := (Y, P) -> not reflq(Y, P);



pnormalform := proc(L, P)
local X, Y, eps, i, Z, W, a;
    X := L[1];
    Y := L[2];
    eps := L[3];
    Z := [supsort(X), supsort(Y), eps*signature(insort(Y))];
    for i to nops(P) do
        a := (M, i, P) -> map((k, i, P) -> P[i][k], M, i, P);
        W := [supsort(map(a, X, i, P)), supsort(map(a, Y, i, P)),
            eps*signature(insort(map(a, Y, i, P)))];
        if lex(W, Z) then Z := W fi
    od;
    Z
end;

    lexico := (L1, L2) -> L1[1] < L2[1] or L1[1] = L2[1] and L1[2] <= L2[2];


                                                                        
bdryplus := (L, k, P) -> pnormalform([L[1], subsop(k = NULL, L[2]), 
             (-1)^k ], P);

        subsetq := (T, S) -> evalb({op(T)} intersect {op(S)} = {op(S)});

signature := proc(G)
local i, j, x;
    x := 1;
    for i to nops(G) do for j from i + 1 to nops(G) do
            if not lexico(G[i], G[j]) then x := -x fi
        od
    od;
    x
end;

zerolist := proc(n)
local i, ans;
    ans := []; for i to n do ans := [op(ans), 0] od; ans
end;

numbholes := proc(M)
# Find the rank of the first homology group of a graph
local i,numb,answ;
answ := [];
for i from 1 to nops(M) do
    numb := 1-nops({op(map(op,M[i][1]))})+nops(M[i][1]);
    answ := [op(answ),numb]
od;
answ
end;

lowranker := proc(kern, cubesused, imfin, kerfin, cubelist, outfilename)
# Uses data for stable homology of holomorphs to determine
# homology for holomorphs in lower dimensions.
local indmin, indmax, Kern, Kernmod, cubesusedholes, cubelistholes,
      kernnullsp,kernbasis,i,j,k, Imfin, Kerfin, Kerfinmod, kerfinnullsp,imfinbasis,
      kerfinbasis;
      cubesusedholes := numbholes(cubesused);
      cubelistholes := numbholes(cubelist);
      indmax := max(op(cubelistholes))-1;
      indmin := min(op(cubesusedholes));
      Kern := linalg[matrix](kern);
      Imfin := linalg[matrix](imfin);
      Kerfin := linalg[matrix](kerfin);
      writeto(outfilename);
      lprint(`Begin low rank calculations.`);
      writeto(terminal);
      for i from indmax by -1 to indmin do
	appendto(outfilename);
        lprint(`---------`);
        lprint(`Rank is `,i,`.`);
        writeto(terminal);
        Kernmod := Kern;
        for j from nops(kern[1]) by -1 to 1 do
	    if cubesusedholes[j] <= i then
	       if linalg[coldim](Kernmod)=1 then
                  Kernmod := {}
               else
	          Kernmod := linalg[delcols](Kernmod,j..j)
              fi
            fi
        od;
        if Kernmod = {} then
	   kernnullsp := {}
        else
           kernnullsp := linalg[nullspace](transpose(Kernmod))
        fi;
        if kernnullsp = {} then
           appendto(outfilename);
	   lprint(`Homology is zero.`);
           writeto(terminal)
        else
	    kernnullsp := linalg[matrix]([op(kernnullsp)]);
            kernbasis := map(M -> convert(M,list), 
                 linalg[rowspace](linalg[multiply](kernnullsp,Kern)));
	    appendto(outfilename);
            lprint(`Kernbasis is `, kernbasis, `.`);
            writeto(terminal);
            Kerfinmod := Kerfin;
	    for j from nops(kerfin[1]) by -1 to 1 do
		if cubelistholes[j] <= i then
                   if linalg[coldim](Kerfinmod)=1 then
                      Kerfinmod := {}
                   else
		      Kerfinmod := linalg[delcols](Kerfinmod,j..j)
                  fi
		fi
            od;
            if Kerfinmod = {} then
	       kerfinnullsp := {}
            else
	       kerfinnullsp := linalg[nullspace](transpose(Kerfinmod))
            fi;
	    if kerfinnullsp = {} then
               appendto(outfilename);
	       lprint(`No elements in image.`);
               writeto(terminal)
	    else
		kerfinnullsp := linalg[matrix]([op(kerfinnullsp)]);
		imfinbasis := map(M -> convert(M,list),
                    linalg[rowspace](linalg[multiply](kerfinnullsp, Imfin)));
                kerfinbasis := map(M -> convert(M,list),
                    linalg[rowspace](linalg[multiply](kerfinnullsp, Kerfin)));
		appendto(outfilename);
		lprint(`imfinbasis is `, imfinbasis, `.`);
		lprint(`kerfinbasis is `, kerfinbasis, `.`);
		writeto(terminal);
            fi
        fi
      od
end;