##############################################################################
##
## SupportOfMorphism(m) returns the list of positions at which the 
## list m is defined.
##
##############################################################################
SupportOfMorphism:=function(m)
    local s, i;
    s:=[];
    for i in [1..Length(m)] do
        if IsBound(m[i]) then Add(s,i);
        fi;
    od;
    return(s);
end;

##############################################################################
##
## IdentityMorphism(l) returns the identity morphism on the set l.
## A morphism is stored as a list of the images of its values on a set
## of numbers, which form its domain of definition, and are taken to be
## an object in a concrete category. At elements of other objects the
## morphism will be undefined.
##
##############################################################################
IdentityMorphism:=function(l)
    local m, i;
    m:=[];
    for i in l do
        m[i]:=i;
    od;
    return(m);
end;

##############################################################################
##
## Composition(f,g) returns the composition of the functions f and g, expressed
## as lists of their values.
##
##############################################################################
Composition:=function(f,g)
    local i, h;
    h:=[];
    for i in [1..Length(f)] do
        if IsBound(f[i]) then
            if IsBound(g[f[i]]) then
                h[i]:=g[f[i]];
            else return(false);
            fi;
        fi;
    od;
    return(h);
end;

##############################################################################
##
## IsComposable(f,g) returns true if the functions f and g, expressed
## as lists of their values, can be composed and false otherwise.
##
##############################################################################
IsComposable:=function(f,g)
    local i;
    for i in [1..Length(f)] do
        if IsBound(f[i]) then
            if not IsBound(g[f[i]]) then return(false);
            fi;
        fi;
    od;
    return(true);
end;

##############################################################################
##
## Objects(cat) returns the objects of the concrete category cat, as a list
## of sets. At the moment it will not work unless for every object there
## is at least one generator morphism whose support is that object.
##
##############################################################################
Objects:=function(cat)
    local m;
    if IsBound(cat.objects) then return(cat.objects);
    fi;
    cat.objects:=[];
    for m in cat.generators do
        Add(cat.objects,SupportOfMorphism(m));
    od;
    cat.objects:=SSortedList(cat.objects);
    return(cat.objects);
end;

##############################################################################
##
## Origin(cat,m) returns the position in cat.objects of the domain of the
## morphism m.
##
##############################################################################
Origin:=function(cat,m)
    local k, x;
    k:=PositionBound(m);
    for x in [1..Length(cat.objects)] do
         if k in cat.objects[x] then
              return(x);
         fi;
    od;
end;

##############################################################################
##
## Terminus(cat,m) returns the position in cat.objects of the domain of the
## morphism m.
##
##############################################################################
Terminus:=function(cat,m)
    local k, x;
    k:=m[PositionBound(m)];
    for x in [1..Length(cat.objects)] do
         if k in cat.objects[x] then
              return(x);
         fi;
    od;
end;



Domains:=function(C)
Objects(C);
C.domain:=List(C.generators,x->Origin(C,x));
C.codomain:=List(C.generators,x->Terminus(C,x));
end;

##############################################################################

ConcreteCategoryOps:=rec();

##############################################################################
##
## ConcreteCategory(list of functions, (list of sets))
## There are optionally one or two arguments. The first is a list of generating
## functions of the category, the second is a list of the objects of the category.
## The function starts a record for a concrete category.
## If there is only one argument, the objects are taken to be the domains of the
## generator morphisms, so for every object there should be
## at least one generator morphism whose support is that object.
## It could be the identity morphisms, but doesn't have to be.
##
## Written by Peter Webb 2008, Moriah Elkin 2018.
##
##############################################################################

ConcreteCategory:=function(args)
    local output, domains, codomains, x, m, nums;
    output:=rec(generators:=args[1], operations:=ConcreteCategoryOps);
    if Length(args)=2 then
        #Checks if object sets overlap or entries repeated (catches ([1,2],[2,3,4]) or ([1,2,2]))
        nums:=[];
        for x in args[2] do
            for m in x do
                if not (m in nums) then
                    Add(nums,m);
                else
                    Error("One or more entries is duplicated within one or more objects.");
                fi;
            od;
        od;

        #Computes domains/codomains of morphisms
        domains:=[];
        for m in args[1] do
            Add(domains,SupportOfMorphism(m));
        od;
        domains:=SSortedList(domains);
        codomains:=[];
        for x in domains do
            for m in output.generators do
                if SupportOfMorphism(m)=x then
                    Add(codomains, List(x,a->m[a]));
                fi;
            od;
        od;
        codomains:=SSortedList(codomains);
        
        #Ensures domains/codomains of morphisms in provided objects list
        if not IsSubset(args[2],domains) then
            Error("One or more morphisms have domains not included in objects list.");
        fi;
        if not IsSubset(args[2],codomains) then
            Error("One or more morphisms have codomains not included in objects list.");
        fi;

        output.objects:=SSortedList(args[2]);
    elif Length(args)=1 then
        Objects(output);
    fi;
    Domains(output);
    return(output);
end;


##############################################################################
##
## EmptyMat(r,s) returns an r x s matrix in which each entry is [].
##
##############################################################################
EmptyMat:=function(r,s)
    local mat, i, j;
    mat:=[];
    for i in [1..r] do
    mat[i]:=[];
        for j in [1..s] do
            mat[i][j]:=[];
        od;
    od;
    return(mat);
end;

##############################################################################
##
## Morphisms(cat) returns an l x l matrix, where l is the number of objects
## in the category cat, and where the i,j entry is a list of the
## morphisms from object i to
## object j.
##
##############################################################################
Morphisms:=function(cat)
    local n, genmat, g, mormat, oldlength, newlength, i, j, k, h, templist;
    if IsBound(cat.morphisms) then return(cat.morphisms);
    fi;
    if not IsBound(cat.objects) then Objects(cat);
    fi;
    n:=Length(cat.objects);
    genmat:=EmptyMat(n,n);
    mormat:=EmptyMat(n,n);
    for g in cat.generators do
        Add(genmat[Origin(cat,g)][Terminus(cat,g)],g);
    od;
    for i in [1..n] do
        Add(mormat[i][i],IdentityMorphism(cat.objects[i]));
    od;
    oldlength:=0;
    newlength:=Sum(List(mormat,x->Sum(List(x,Length))));
    while oldlength < newlength do
        oldlength:=newlength;
        for i in [1..n] do
            for j in [1..n] do
                templist:=[];
                for k in [1..n] do
                    for g in genmat[k][j] do
                        for h in mormat[i][k] do
                            Add(templist, Composition(h,g));
                        od;
                    od;
                 od;
                 Append(mormat[i][j],templist);
                 mormat[i][j]:=SSortedList(mormat[i][j]);
            od;
        od;
        newlength:=Sum(List(mormat,x->Sum(List(x,Length))));
    od;
    cat.morphisms:=mormat;
    return(mormat);
end;


##############################################################################
##
## Chains(cat,n) returns an l x l matrix, where l is the number of objects
## in the category cat, and where the i,j entry is a list of the chains
## of length n of composable morphisms, starting at object i and finishing
## at object j. The degenerate chains (i.e. those with an identity morphism
## somewhere in the list) are included.
##
##############################################################################
Chains:=function(cat,n)
local s, i, j, k, l, u, v, newchain;
    if not IsBound(cat.chains) then
        if not IsBound(cat.morphisms) then Morphisms(cat);
        fi;
        cat.chains:=[]; cat.chains[1]:=List(cat.morphisms,x->List(x,y->List(y,z->[z])));
    fi;
    if IsBound(cat.chains[n]) then return(cat.chains[n]);
    fi;
    s:=Length(cat.objects);
    while Length(cat.chains) < n do
        i:=Length(cat.chains);
        cat.chains[i+1]:=EmptyMat(s,s);
        for j in [1..s] do
            for k in [1..s] do
                for l in [1..s] do
                    for u in cat.chains[i][j][l] do
                        for v in cat.morphisms[l][k] do
                            newchain:=[];
                            Append(newchain,u);
                            Add(newchain,v);
                            Add(cat.chains[i+1][j][k],newchain);
                        od;
                    od;
                od;
                cat.chains[i+1][j][k]:=SSortedList(cat.chains[i+1][j][k]);

            od;
        od;
    od;
return(cat.chains[n]);
end;


##############################################################################
##
## NDChains(cat,n) returns an l x l matrix, where l is the number of objects
## in the category cat, and where the i,j entry is a list of the non-degenerate
## chains of length n of composable morphisms, starting at object i and finishing
## at object j.
##
##############################################################################
NDChains:=function(cat,n)
local s, i, j, k, l, m, u, v, newchain, position;
    if not IsBound(cat.morphisms) then Morphisms(cat);
    fi;
    s:=Length(cat.objects);
    if not IsBound(cat.ndchains) then
        cat.ndchains:=[];
        cat.ndchainpositions:=[];
        cat.ndchains[1]:=EmptyMat(s,s);
        cat.ndchainpositions[1]:=EmptyMat(s,s);
        position:=1;
        for j in [1..s] do
            for k in [1..s] do
                for v in cat.morphisms[j][k] do
                    if j=k then
                        m:=IdentityMorphism(cat.objects[j]);
                        if v<>m then Add(cat.ndchains[1][j][k],[v]);
                        Add(cat.ndchainpositions[1][j][k],position);
                        position:=position+1;
                        fi;
                    else Add(cat.ndchains[1][j][k],[v]);
                    Add(cat.ndchainpositions[1][j][k],position);
                    position:=position+1;
                    fi;
                 od;
            od;
        od;
    fi;
    if IsBound(cat.ndchains[n]) then return(cat.ndchains[n]);
    fi;
    while Length(cat.ndchains) < n do
        i:=Length(cat.ndchains);
        cat.ndchains[i+1]:=EmptyMat(s,s);
        cat.ndchainpositions[i+1]:=EmptyMat(s,s);
        position:=1;
        for j in [1..s] do
            for k in [1..s] do
                for l in [1..s] do
                    for u in cat.ndchains[i][j][l] do
                        for v in cat.ndchains[1][l][k] do
                            newchain:=[];
                            Append(newchain,u);
                            Append(newchain,v);
                            Add(cat.ndchains[i+1][j][k],newchain);
                            Add(cat.ndchainpositions[i+1][j][k],position);
                            position:=position+1;
                        od;
                    od;
                od;
                cat.ndchains[i+1][j][k]:=SSortedList(cat.ndchains[i+1][j][k]);

            od;
        od;
    od;
return(cat.ndchains[n]);
end;



##############################################################################
##
## BoundaryComponent(cat,n,i) returns the integer valued matrix which represents
## component i of the boundary map Zcat.chains[n] -> Zcat.chains[n-1]. The
## integer i must lie between 0 and n.
## This component sends a chain [m1,m2, ... ,mn]
## to [m2, ... ,mn] if i=0,
## to [m1, ... , mi*mi+1, ... ,mn] if i lies between 1 and n-1, and
## to [m1, ... , mn-1] if i=n.
## All chains are used, including the degenerate ones.
##
##############################################################################
BoundaryComponent:=function(cat,n,i)
    local s, matrix, image, u, v, w, z, newchain, p, q;
    Chains(cat,n);
    s:=Length(cat.objects);
    matrix:=[];
    if n=1 then
        if i=0 then
            image:=List(cat.objects,x->0);
            for u in [1..s] do
                for v in [1..s] do
                    for z in cat.morphisms[u][v] do
                        image[v]:=1;
                        Add(matrix,ShallowCopy(image));
                        image[v]:=0;
                    od;
                od;
            od;
        fi;
        if i=1 then
            image:=List(cat.objects,x->0);
            for u in [1..s] do
                for v in [1..s] do
                    for z in cat.morphisms[u][v] do
                        image[u]:=1;
                        Add(matrix,ShallowCopy(image));
                        image[u]:=0;
                    od;
                od;
            od;
        fi;
        return(matrix);
    fi;
    image:=List(cat.chains[n-1],x->List(x,y->List(y,w->0)));
    for u in [1..s] do
        for v in [1..s] do
            for z in cat.chains[n][u][v] do
                newchain:=[];
                if i=0 then
                    newchain:=z{[2..n]};
                    q:=Origin(cat,newchain[1]);
                    p:=Position(cat.chains[n-1][q][v],newchain);
                    image[q][v][p]:=1;
                    Add(matrix,Flat(image));
                    image[q][v][p]:=0;
                elif i=n then
                    newchain:=z{[1..n-1]};
                    q:=Terminus(cat,newchain[n-1]);
                    p:=Position(cat.chains[n-1][u][q],newchain);
                    image[u][q][p]:=1;
                    Add(matrix,Flat(image));
                    image[u][q][p]:=0;
                else
                    newchain:=z{[1..i-1]};
                    newchain[i]:=Composition(z[i],z[i+1]);
                    Append(newchain,z{[i+2..n]});
                    p:=Position(cat.chains[n-1][u][v],newchain);
                    image[u][v][p]:=1;
                    Add(matrix,Flat(image));
                    image[u][v][p]:=0;
                fi;
            od;
        od;
    od;
    return(matrix);
end;

##############################################################################
##
## NDBoundaryComponent(cat,n,i) returns the integer valued matrix which represents
## component i of the boundary map Zcat.ndchains[n] -> Zcat.ndchains[n-1]. The
## integer i must lie between 0 and n.
## This component sends a chain [m1,m2, ... ,mn]
## to [m2, ... ,mn] if i=0,
## to [m1, ... , mi*mi+1, ... ,mn] if i lies between 1 and n-1, and
## to [m1, ... , mn-1] if i=n
## where these symbols are zero if any morphism is the identity.
## Only non-degenerate chains are used.
##
##############################################################################
NDBoundaryComponent:=function(cat,n,i)
    local s, matrix, image, u, v, w, z, newchain, p, q;
    NDChains(cat,n);
    s:=Length(cat.objects);
    matrix:=[];
    if n=1 then
        if i=0 then
            image:=List(cat.objects,x->0);
            for u in [1..s] do
                for v in [1..s] do
                    for z in cat.ndchains[1][u][v] do
                        image[v]:=1;
                        Add(matrix,ShallowCopy(image));
                        image[v]:=0;
                    od;
                od;
            od;
        fi;
        if i=1 then
            image:=List(cat.objects,x->0);
            for u in [1..s] do
                for v in [1..s] do
                    for z in cat.ndchains[1][u][v] do
                        image[u]:=1;
                        Add(matrix,ShallowCopy(image));
                        image[u]:=0;
                    od;
                od;
            od;
        fi;
        return(matrix);
    fi;
    image:=List(cat.ndchains[n-1],x->List(x,y->List(y,w->0)));
    for u in [1..s] do
        for v in [1..s] do
            for z in cat.ndchains[n][u][v] do
                newchain:=[];
                if i=0 then
                    newchain:=z{[2..n]};
                    q:=Origin(cat,newchain[1]);
                    p:=Position(cat.ndchains[n-1][q][v],newchain);
                    image[q][v][p]:=1;
                    Add(matrix,Flat(image));
                    image[q][v][p]:=0;
                elif i=n then
                    newchain:=z{[1..n-1]};
                    q:=Terminus(cat,newchain[n-1]);
                    p:=Position(cat.ndchains[n-1][u][q],newchain);
                    image[u][q][p]:=1;
                    Add(matrix,Flat(image));
                    image[u][q][p]:=0;
                else
                    newchain:=z{[1..i-1]};
                    newchain[i]:=Composition(z[i],z[i+1]);
                    if newchain[i]=IdentityMorphism(cat.objects[Origin(cat,newchain[i])]) then
                        Add(matrix,Flat(image));
                    else
                        Append(newchain,z{[i+2..n]});
                        p:=Position(cat.ndchains[n-1][u][v],newchain);
                        image[u][v][p]:=1;
                        Add(matrix,Flat(image));
                        image[u][v][p]:=0;
                    fi;
                fi;
            od;
        od;
    od;
    return(matrix);
end;


##############################################################################
##
## BoundaryMap(cat,n) returns the integer valued matrix which represents
## the boundary map Zcat.chains[n] -> Zcat.chains[n-1].
## It is Sum (-1)^i BoundaryComponent(cat,n,i)
## All chains are used, including the degenerate ones.
##
##############################################################################
BoundaryMap:=function(cat,n)
    local mat, i;
    mat:=BoundaryComponent(cat,n,0);
    for i in [1..n] do
        mat:= mat + (-1)^i*BoundaryComponent(cat,n,i);
    od;
    return(mat);
end;

##############################################################################
##
## NDBoundaryMap(cat,n) returns the integer valued matrix which represents
## the boundary map Zcat.ndchains[n] -> Zcat.ndchains[n-1].
## It is Sum (-1)^i NDBoundaryComponent(cat,n,i)
## Only non-degenerate chains are used.
##
##############################################################################
NDBoundaryMap:=function(cat,n)
    local mat, i;
    mat:=NDBoundaryComponent(cat,n,0);
    for i in [1..n] do
        mat:= mat + (-1)^i*NDBoundaryComponent(cat,n,i);
    od;
    return(mat);
end;

##############################################################################
##
## NDChainImage(cat,c) returns the integer valued vector which represents
## image of the chain c under the boundary map
## Zcat.ndchains[n] -> Zcat.ndchains[n-1] where n is the length of c. 
## The map sends a chain c=[m1,m2, ... ,mn]
## to [m2, ... ,mn] + Sum (-1)^i*[m1, ... , mi*mi+1, ... ,mn]
## + (-1)^n*[m1, ... , mn-1]
## where these symbols are zero if any morphism is the identity.
## Only non-degenerate chains are used.
##
##############################################################################
NDChainImage:=function(cat,c)
local n, u, v, image, newchain, q, p, vec, i;
    n:=Length(c);u:=Origin(cat,c[1]);v:=Terminus(cat,c[n]);
    if n<2 then Error("The chain must have length at least 2.");
    fi;
    image:=List(cat.ndchains[n-1],x->List(x,y->List(y,w->0)));
    newchain:=[];
    newchain:=c{[2..n]};
    q:=Origin(cat,newchain[1]);
    p:=Position(cat.ndchains[n-1][q][v],newchain);
    image[q][v][p]:=1;
    vec:=Flat(image);
    image[q][v][p]:=0;
    for i in [1..n-1] do
         newchain:=c{[1..i-1]};
         newchain[i]:=Composition(c[i],c[i+1]);
         if newchain[i]<>IdentityMorphism(cat.objects[Origin(cat,newchain[i])]) then
              Append(newchain,c{[i+2..n]});
              p:=Position(cat.ndchains[n-1][u][v],newchain);
              image[u][v][p]:=1;
              vec:=vec+(-1)^i*Flat(image);
              image[u][v][p]:=0;
         fi;
    od;
    newchain:=c{[1..n-1]};
    q:=Terminus(cat,newchain[n-1]);
    p:=Position(cat.ndchains[n-1][u][q],newchain);
    image[u][q][p]:=1;
    vec:=vec+(-1)^n*Flat(image);
    return(vec);
end;

##############################################################################
##
## NDShortChainImage(cat,c) returns the image of the chain c under the
## boundary map, stored in a short form.
## Writing c=[m1,m2, ... ,mn] the short form is the list of the positions
## of the n+1 chains
## [m2, ... ,mn], ..., [m1, ... , mi*mi+1, ... ,mn], ..., [m1, ... , mn-1]
## in cat.ndchains[n-1]. The position is the corresponding number in
## cat.ndchainpositions[n-1].
## If any of the products computed here is the identity, the corresponding
## entry is left unbound. The length of c must be at least 2.
## Only non-degenerate chains are used.
## The idea behind the short form is to exploit the sparseness of the matrix
## of the boundary map, thus saving space. Because the code is not compiled,
## there may only be a saving of time compared to routines which use
## built-in matrix operations of GAP when the matrix is large.
##
##############################################################################


NDShortChainImage:=function(cat,c)
local n, u, v, image, newchain, q, i, position;
    n:=Length(c);u:=Origin(cat,c[1]);v:=Terminus(cat,c[n]);
    if n<2 then Error("The chain must have length at least 2.");
    fi;
    newchain:=[];
    image:=[];
    newchain:=c{[2..n]};
    q:=Origin(cat,newchain[1]);
    position:=Position(cat.ndchains[n-1][q][v],newchain);
    image[1]:=cat.ndchainpositions[n-1][q][v][position];
    for i in [1..n-1] do
         newchain:=c{[1..i-1]};
         newchain[i]:=Composition(c[i],c[i+1]);
         if newchain[i]<>IdentityMorphism(cat.objects[Origin(cat,newchain[i])])
              then Append(newchain,c{[i+2..n]});
              position:=Position(cat.ndchains[n-1][u][v],newchain);
              image[i+1]:=cat.ndchainpositions[n-1][u][v][position];
         fi;
    od;
    newchain:=c{[1..n-1]};
    q:=Terminus(cat,newchain[n-1]);
    position:=Position(cat.ndchains[n-1][u][q],newchain);
    image[n+1]:=cat.ndchainpositions[n-1][u][q][position];
    return(image);
end;

##############################################################################
##
## NDShortCheckRoutine(cat,c,check,p) takes a non-degenerate n-chain c
## and a matrix 'check' whose
## number of columns is the number of (n-1)-chains, and returns a matrix whose
## rows span the subspace of the span of 'check' orthogonal to the boundary
## d(c). If the rows of check are independent, then so are the rows of the
## output.
## The algorithm exploits the sparseness of d(c) (it has n+1 non-zero entries),
## and also achieves speed by only operating on the rows of check which are
## not already orthogonal to d(c). This check of orthogonality is only a few
## multiplications and is quick.
## The algorithm works over GF(p). The chain c must have length at least 2.
##
##############################################################################

NDShortCheckRoutine:=function(cat,c,check,p)
    local firstfound, firstposition, l, image, v, product, i, pivotvec, temp;
    firstfound:=false;
    if check=[] then return(check);
    fi;
    image:=NDShortChainImage(cat,c);
    if image=[] then return(check);
    fi;
    l:=Length(check);
    for v in [1..l] do
        product:=Zero(GF(p));
        for i in [1..Length(c)+1] do
            if IsBound(image[i]) then
                product:=product+(-1)^(i+1)*check[v][image[i]];
            fi;
        od;
        if product<>Zero(GF(p)) then
            if not firstfound then
                firstfound:=true;
                firstposition:=v;
                pivotvec:=product^(-1)*check[v];
            else
                check[v]:=check[v]-product*pivotvec;
            fi;
        fi;
    od;
    if IsBound(firstposition) then
        check{[firstposition..l-1]}:=check{[firstposition+1..l]};
        Unbind(check[l]);
    fi;
    return(check);
end;

##############################################################################
## 
## NDShortBoundaryCheckMatrix(cat,n,p) returns a matrix over GF(p) whose rows
## are a basis for the orthogonal space to the image of the boundary map
## from non-degenerate n-chains to (n-1)-chains. n must be at least 2.
##
##############################################################################

NDShortBoundaryCheckMatrix:=function(cat,n,p)
local s, t, check, u, v, c;
    NDChains(cat,n);
    s:=Length(cat.objects);
    t:=Sum(Sum(List(cat.ndchains[n-1],x->List(x,Length))));
    check:=IdentityMat(t,GF(p));
    for u in [1..s] do
        for v in [1..s] do
            for c in cat.ndchains[n][u][v] do
                check:=NDShortCheckRoutine(cat,c,check,p);
           od;
        od;
    od;
    return(check);
end;

##############################################################################
##
## Homology(cat,n) returns the invariant factors of the homology of
## the nerve of cat in dimension n.
## Non-degenerate chains are used. The answer [] means the homology is zero.
## Zeros in the answer correspond to infinite cyclic factors.
##
##############################################################################

Homology:=function(cat,n)
    local mat1, l, mat2, factors, i;
    Objects(cat);
    if n=0 then
         l:=Length(cat.objects)-RankMat(NDBoundaryMap(cat,1));
         return(NullMat(1,l)[1]);
    fi;
    mat1:=NDBoundaryMap(cat,n);
    if mat1=[] then return([]);
    fi;
    l:=Length(NullspaceMat(mat1));
    mat2:=NDBoundaryMap(cat,n+1);
    if mat2=[] then return(NullMat(1,l)[1]);
    fi;
##
## The next lines are to deal with an apparent bug in GAP, that
## SmithNormalFormIntegerMat of a zero matrix produces an error.
##
    factors:=ShallowCopy(mat2[1]);
    if mat2 <> NullMat(Length(mat2),Length(mat2[1])) then
        factors:=Sum(SmithNormalFormIntegerMat(NDBoundaryMap(cat,n+1)));
    fi;
    i:=1;
    while factors[i]=1 do
        i:=i+1;
    od;
    if i>l then
        return([]);
    fi;
    return(factors{[i..l]});
end;

##############################################################################
##
## HomologyDimension(cat,n,p) returns the dimension of the homology of
## the nerve of cat in dimension n over GF(p) by using non-degenerate chains.
##
##############################################################################
HomologyDimension:=function(cat,n,p)
    local mat1, mat2;
    Objects(cat);
    if n=0 then
         return(Length(cat.objects)-RankMat(Identity(GF(p))*NDBoundaryMap(cat,1)));
    fi;
    mat1:=Identity(GF(p))*NDBoundaryMap(cat,n);
    if mat1=[] then return(0);
    fi;
    mat2:=Identity(GF(p))*NDBoundaryMap(cat,n+1);
    if mat2=[] then return(Length(mat1)-RankMat(mat1));
    fi;
    return(Length(mat1)-RankMat(mat1)-RankMat(mat2));
end;

##############################################################################
##
## NDShortHomologyDimension(cat,n,p) returns the dimension of the homology of
## the nerve of cat in dimension n over GF(p) by using non-degenerate chains.
## The 'short' space-saving algorithm is used.
##
##############################################################################
NDShortHomologyDimension:=function(cat,n,p)
    local a, b, c;
    Objects(cat);
    if n=0 then
         return(Length(cat.objects)-RankMat(Identity(GF(p))*NDBoundaryMap(cat,1)));
    fi;
    NDChains(cat,n);
    if n=1 then
        return(Length(NDShortBoundaryCheckMatrix(cat,n+1,p))
        -RankMat(Identity(GF(p))*NDBoundaryMap(cat,1)));
    fi;
     if n>1 then
        a:=Sum(Sum(List(cat.ndchains[n-1],x->List(x,Length))));
        b:=Length(NDShortBoundaryCheckMatrix(cat,n,p));
        c:=Length(NDShortBoundaryCheckMatrix(cat,n+1,p));
        return(b+c-a);
    fi;
end;
##############################################################################
##
## The next two routines produce an initial check matrix consisting of a basis
## for the orthogonal space to the images of a set of random chains. The idea
## is to cut down the size of the check matrix as soon as possible. This only
## need be done with very large matrices.
##
##############################################################################

RandomChains:=function(cat,n)
    local s, chains, i, j, l, m, indices, r;
    NDChains(cat,n);
    s:=Length(cat.objects);
    chains:=[];
    for i in [1..s] do
        for j in [1..s] do
            l:=Length(cat.ndchains[n-1][i][j]);
            m:=Length(cat.ndchains[n][i][j]);
            indices:=[];
            while Length(indices)<l do
                r:=Random([1..m]);
                if not r in indices then Add(indices,r);
                fi;
            od;
            Append(chains,cat.ndchains[n][i][j]{Set(indices)});
        od;
    od;
    return(chains);
end;
StarterCheckMatrix:=function(cat,chains,p)
local n, t, check, c;
    n:=Length(chains[1]);
    t:=Sum(Sum(List(cat.ndchains[n-1],x->List(x,Length))));
    check:=IdentityMat(t,GF(p));
    for c in chains do
        check:=NDShortCheckRoutine(cat,c,check,p);
    od;
    return(check);
end;