Quelle  subdivision.gi   Sprache: unbekannt

#######################################################################
##
##  RigidFacetsSubdivision
##  Input: A non-rigid cell complex C, with finite cell stabilizers and such
##  that we could equip it with a metric such that each of the cells of C is
##  convex, the restriction of the metric to each cell is CAT(0)
##  and the cell stabilizers act by CAT(0) isometries.
##  An optional second argument allows to carry out the subdivision
##  only up to a given upper bound on the cell dimension.
##
##  Output: A rigidification of C.
##
####################################################################
DeclareGlobalFunction("ListsOfCellsToRegularCWComplex");
InstallGlobalFunction(ListsOfCellsToRegularCWComplex,
function(dim,boundaries,coboundaries,orientation)
local NrCells,Properties;

Properties:=[["dimension",dim]];

NrCells:=function(n);
if n>dim then return 0; fi;
return Length(Filtered(boundaries[n+1],x->not x[1]=0));
end;

return Objectify(HapRegularCWComplex,
       rec(
           nrCells:=NrCells,
           boundaries:=boundaries,
           coboundaries:=coboundaries,
           orientation:=orientation,
           vectorField:=fail,
           inverseVectorField:=fail,
           criticalCells:=fail,
           properties:=Properties));

end);
##################################################################
InstallGlobalFunction(RigidFacetsSubdivision,
function(arg)
local W, StabRec, i, j, N, M, x, bdry, s1, s2, p, k, w, t,
      DimRec, BoundaryRec, id, dims, NotRigid, NewCell,
      Cell, Elts, Boundary, Dimension, CLeftCosetElt, LCoset,
      pos, Stab, Mult,DimTemp, BoundaryTemp, Partition, ConnectedCheck,
      Stabilizer, Action, IsRigidCell, ReplaceCell, SubdividingCell,
      Orbit, C, NaiveEulerChar, CellHomology, exportFile, IsContractible,
      CosetRec, IsSimplyConnected, BoundaryMult, BoundaryMultRec, jj, jjj;

    C:=arg[1];
    i:=0;
    while C!.dimension(i)>0 do
        i:=i+1;
    od;
    N:=i-1; # Length of the chain complex
    M:=N;


if Length(arg)=2 then M:=Minimum(N,arg[2]);
fi;


    Elts:=C!.elts;
    StabRec:=[];
    DimRec:=[];
    CosetRec:=[];

    ##################################################################
    # If g in Elts return the position of g in the list,
    # otherwise, add g to Elts and return the position.
    pos:=function(g)
    local posit;

        posit:=Position(Elts,g);
        if posit=fail then
            Add(Elts,g);
            return Length(Elts);
        else
            return posit;
        fi;
    end;
    ##################################################################
    id:=pos(One(C!.group));
    ##################################################################
    # return the stabilizer of g*e,
    #
    Stab:=function(e,g)
    return ConjugateGroup(StabRec[e[1]+1][AbsInt(e[2])],Elts[g]^-1);
    end;
    ##################################################################
    # returns  a  "canonical"  representative  of  the  right  coset
    # Elts[g]*Stab[i+1][j]

    ######Add a record for left coset
    for i in [0..N] do
      CosetRec[i+1]:=[];
    od;

    CLeftCosetElt:=function(i,j,g)
      local p;
      if IsBound(CosetRec[AbsInt(i)+1][j]) then
        if IsBound(CosetRec[AbsInt(i)+1][j][g]) then
          return CosetRec[AbsInt(i)+1][j][g];
        else
          p:=pos(CanonicalRightCountableCosetElement
                                  (StabRec[AbsInt(i)+1][j],Elts[g]^-1)^-1);
          CosetRec[AbsInt(i)+1][j][g]:=p;
          return p;
        fi;
      else
        CosetRec[i+1][j]:=[];
        p:=pos(CanonicalRightCountableCosetElement
                                (StabRec[AbsInt(i)+1][j],Elts[g]^-1)^-1);
        CosetRec[AbsInt(i)+1][j][g]:=p;
        return p;
      fi;

    end;
    ##################################################################
    ##
    ##  Input:  A list L, degree k, position g of an element
    ##  Output: Product of g and L.
    ##
    Mult:=function(L,k,g)
    local x,w,t,h,y,vv,LL;
        vv:=[];
        LL:=ShallowCopy(L);
        for x in [1..Length(LL)] do
            w:=Elts[g]*Elts[LL[x][2]];
            t:=CLeftCosetElt(k,AbsInt(LL[x][1]),pos(w));
            Add(vv,[LL[x][1],t]);
        od;
        return vv;
    end;
    ###################################################################
    # Store essential data: stabilizers, boundaries, dimensions

    Orbit:=[];
    for i in [1..N+1] do
        Orbit[i]:=[];
    od;

    for i in [0..N] do
        StabRec[i+1]:=[];
        DimRec[i+1]:=C!.dimension(i);
        for j in [1..C!.dimension(i)] do
            StabRec[i+1][j]:=C!.stabilizer(i,j);
        od;
    od;

    BoundaryRec:=[];
    for i in [1..N] do
        BoundaryRec[i]:=[];
        for j in [1..DimRec[i+1]] do
            bdry:=C!.boundary(i,j);
            BoundaryRec[i][j]:=[];
            for x in bdry do
                s1:=C!.action(i-1,AbsInt(x[1]),x[2]);
                p:=pos(CanonicalRightCountableCosetElement
                            (C!.stabilizer(i-1,AbsInt(x[1])),Elts[x[2]]^-1)^-1);
                s2:=C!.action(i-1,AbsInt(x[1]),p);

                Add(BoundaryRec[i][j],[s1*s2*x[1],p]);
            od;

        od;
    od;

    #############
    BoundaryMultRec:=[];
    for jj in [1..N] do
      BoundaryMultRec[jj]:=[];
      for jjj in [1..Length(BoundaryRec[jj])] do
        BoundaryMultRec[jj][jjj]:=[];
      od;
    od;
    #############
    BoundaryMult:=function(i,j,g)
      local w;
      if IsBound(BoundaryMultRec[i]) then
        if IsBound(BoundaryMultRec[i][j]) then
          if IsBound(BoundaryMultRec[i][j][g]) then
            return BoundaryMultRec[i][j][g];
          else
            BoundaryMultRec[i][j][g]:=Mult(BoundaryRec[i][j],i-1,g);
            return BoundaryMultRec[i][j][g];
          fi;
        else
          BoundaryMultRec[i][j]:=[];
          BoundaryMultRec[i][j][g]:=Mult(BoundaryRec[i][j],i-1,g);
          return BoundaryMultRec[i][j][g];
        fi;
      else
        BoundaryMultRec[i]:=[];
        BoundaryMultRec[i][j]:=[];
        BoundaryMultRec[i][j][g]:=Mult(BoundaryRec[i][j],g);
        return BoundaryMultRec[i][j][g];
      fi;
    end;
    ##################################################################


    ##################################################################
    # Check whether the cell is rigid or not

    IsRigidCell:=function(k,m)
    local bdry, intst, L;
        bdry:=BoundaryRec[k][m];
        L:=List(bdry,w->Elements(ConjugateGroup(StabRec[k][AbsInt(w[1])],Elts[w[2]]^-1)));
        intst:=Intersection(L);
        if not Elements(StabRec[k+1][m])=Elements(intst) then
            return false;
        else return true;
        fi;

    end;
    ##################################################################


    ##################################################################
    # Subdividing (k+1)-cell number i (denoted by sigma) with the
    # Rigid Facets Subdivision Algorithm :
    SubdividingCell:=function(k,i)
    local bdry, w, x, d, y, a, b, T, j, j1,j2, z, t, c, Flag, s, s1, ConnectToCenter,
          SearchComponent,temp, l1, bdry1, bdry2, Orb, IsSameOrbit, OrbFlag, OrbElm,
          ConnectedCheck,tau,Ctau,F,exportBoundary,exportStabilizer,
          Boundaries,Stabilizers,Tcells, convexTcells,p;

 # Record the barycenter of sigma as a vertex, with stabilizer Gamma_sigma :
        DimRec[1]:=DimRec[1]+1;
        Add(StabRec[1],StabRec[k+1][i]);

        bdry:=ShallowCopy(BoundaryRec[k][i]);


        #############################################################
        # In each orbit {g_{jt} sigma_j with j >= 2,
 # choose one cell such that its union with the cells in T is connected,
 # and such that the naive Euler characteristic of this union remains 1.
  SearchComponent:=function(q)
  local x,z,j,w,TT, tau, boundaryOFtau, S, Ps, ps,countPrint, bool;
      x:=[];
      w:=[];
      Ps:=[[1,1]];

      j:=1;
      while not Length(T)=Length(Orb) do


          while j<=Length(Orb) do
              if OrbFlag[j]=0 then

                  while j1<=Length(Orb[j]) do
                      if Flag[j][j1]=0 then
                          TT:=Union(T,[Orb[j][j1]]);



                          if IsContractible(k-1,TT) then

                            if k>=2 then
                              S:=[];
                              for tau in TT do

                                   boundaryOFtau:=BoundaryMult(k-1,AbsInt(tau[1]),tau[2]);
                                if  tau[1]<0 then
                                    boundaryOFtau:=NegateWord(boundaryOFtau);
                                fi;
                                Append(S,boundaryOFtau);

                              od;
                              S:=AlgebraicReduction(S);

                              if (Length(T)=Length(Orb)-1) and k>=4 then bool:=true;
                              else bool:=false;
                              fi;
                              if IsSimplyConnected(k-2,S,bool) then

                                  Add(T,Orb[j][j1]);
                                  Flag[j][j1]:=1;
                                  Add(Ps,[j,j1]);
                                  OrbFlag[j]:=1;
                                  j:=1;
                                  break;
                              fi;
                            else
                              Add(T,Orb[j][j1]);
                              Flag[j][j1]:=1;
                              Add(Ps,[j,j1]);
                              OrbFlag[j]:=1;
                              j:=1;
                              break;

                            fi;

                          fi;

                      fi;

                      j1:=j1+1;
                    od;
                  if j1<=Length(Orb[j]) then
                      break;

                  fi;
              fi;
              j:=j+1;
              j1:=1;
          od;
          if j>Length(Orb) then
              Remove(T);

              ps:=Ps[Length(Ps)];
              OrbFlag[ps[1]]:=0;
              Flag[ps[1]][ps[2]]:=0;
              Remove(Ps);
              j:=ps[1];
              j1:=ps[2]+1;
              while j1>Length(Orb[j]) do
                ps:=Ps[Length(Ps)];
                OrbFlag[ps[1]]:=0;
                Flag[ps[1]][ps[2]]:=0;
                Remove(Ps);
                Remove(T);
                j:=ps[1];
                j1:=ps[2]+1;

              od;
              if Length(T)=0 then
                  Print("\n Could not find any fundamental domain for the ",i,"th of the",k,"-cells! \n");
                  Print(1/0);
              fi;
          fi;

      od;

      return T;
  end;
        #############################################################

    #################################################################
    # Connect the collection of (k-1)-cells T to the barycenter of the cell sigma.
    # The elements of T, as well as sigma, are in the format [k,i,gamma]:
    # dimension k, obtained by sending the ith orbit representative to its image under the element gamma (specified by its number in Gamma via C!.elts).
    # Input: e := [k-1, T].
    ConnectToCenter:=function(e)
    local bdry1, x, stab, bdrye, w, stablst, redbdry, tau, boundaryOFtau,
          LCoset, AddCell;


    ##################################################################
    # Add a k-cell with stabilizer stab and boundary bdry
    # to the cell complex
    AddCell:=function(m,stab,bdry,e)
    local a,j,g,s,w,u,v;

    if m<k then
        w:=e[1];

        for j in [1..DimTemp[m+1]] do
            u:=BoundaryTemp[m+1][j];
     # The cells u and w are in the format [orbit number, gamma],
     # such that gamma (pointing to C!.elts) sends the orbit representative to the cell.
     # We make use of the common vertex, namely the barycenter of sigma,
     # to test if there exists gamma in Gamma_sigma with gamma u = w.
            s:=DimRec[m+1]-DimTemp[m+1]+j;
            if AbsInt(w[1])=AbsInt(u[1]) then

             # v:=StabRec[m][AbsInt(u[1]);
              for v in StabRec[m][AbsInt(u[1])] do
                g:=Elts[w[2]]*v*Elts[u[2]]^-1;
                 if g in StabRec[k+1][i] then
                    return [s, CLeftCosetElt(m,s,pos(g))];
                  fi;
              od;

            fi;
        od;

        DimTemp[m+1]:=DimTemp[m+1]+1;

        Add(BoundaryTemp[m+1],e[1]);

 #Update the dimension :
        DimRec[m+1]:=DimRec[m+1]+1;

 # Record the stabilizer :
        Add(StabRec[m+1],stab);\

 # Record the boundary of F :
        Add(BoundaryRec[m],bdry);

    else
        DimRec[m+1]:=DimRec[m+1]+1;
        Add(StabRec[m+1],stab);
         Add(BoundaryRec[m],bdry);
    fi;
    return [DimRec[m+1],CLeftCosetElt(m,DimRec[m+1],id)];
    end;
    ##################################################################

        if e[1]=0 # k-1 = 0 : e[2] = T is a collection of vertices.
 then
            p:=Position(Tcells[1],e[2][1]);
            if not p=fail then return convexTcells[1][p];
            else
            bdry1:=[[-SignInt(e[2][1][1])*DimRec[1],CLeftCosetElt(0,DimRec[1],id)],[e[2][1][1],e[2][1][2]]];

            stab:=Intersection(Stab([0,e[2][1][1]],e[2][1][2]),StabRec[k+1][i]);
            w:=AddCell(1,stab,bdry1,e[2]);
            Add(Tcells[1],e[2][1]);
            Add(convexTcells[1],w);
            return w;
            fi;
        fi;

        stablst:=[];
        bdry1:=[];
        Append(bdry1,e[2]);
        bdrye:=[];

 ## Enumerate the set S := boundary(geometric realization of T) =: bdrye,
 ## by letting tau run through the cells of T = e[2],
 ## and afterwards cancelling cells in S which appear in pairs with opposite orientation.
        for tau in e[2] do

            boundaryOFtau:=BoundaryMult(e[1],AbsInt(tau[1]),tau[2]);
            if  tau[1]<0 then boundaryOFtau:=NegateWord(boundaryOFtau);fi;
            Append(bdrye,boundaryOFtau);
            Add(stablst,Stab([e[1],tau[1]],tau[2]));
        od;
        Add(stablst,StabRec[k+1][i]);
 ## Cancel cells in S which appear in pairs with opposite orientation :
        bdrye:=AlgebraicReduction(bdrye);
 ## Construction of S = bdrye is now completed.

 ## For every x in S, take the convex envelope w := e(x) of x and the barycenter of sigma :

        for x in bdrye do
          p:=Position(Tcells[e[1]],[AbsInt(x[1]),x[2]]);
          if not p=fail then
            w:=convexTcells[e[1]][p];
          else
            w:=ConnectToCenter([e[1]-1,[[AbsInt(x[1]),x[2]]]]);
            Add(Tcells[e[1]],[AbsInt(x[1]),x[2]]);
            Add(convexTcells[e[1]],w);
          fi;
            Add(bdry1,[-SignInt(x[1])*w[1],w[2]]);
     ## Here, the orientation of w := e(x) has been recorded in -SignInt(x[1]).
        od;

 # Calculate the stabilizer :
        stab:=Intersection(stablst);
        if not Length(e[2])=1 then

            w:=AddCell(e[1]+1,stab,bdry1,e[2]);

            return w;
        else
            return AddCell(e[1]+1,stab,bdry1,e[2]);
        fi;
    end;
    ##################################################################
    ## End of function ConnectToCenter, end of Algorithm 3.
    ##################################################################


    IsSameOrbit:=function(e,f)
    local s1,s;
        if AbsInt(e[1])=AbsInt(f[1]) then
            s:=StabRec[k+1][i];
            s1:=StabRec[k][AbsInt(e[1])];
            for x in s do
                if Elts[f[2]]^-1*x*Elts[e[2]] in s1 then
                    return SignInt(e[1])*SignInt(f[1])*pos(x);
                fi;
            od;
            return false;
         else
            return false;
         fi;
    end;

##### 2017 - LUXEMBOURG


##### Check if the fundamental domain formed by sub[1] is connected ####

ConnectedCheck:=function(F,n)
 local b, lst, count, flag, i, d, cnr, m;

 ###

 lst:=[];
        flag:=[];
 b:=[];
 for i in [1..Length(F)] do
  flag[i]:=0;
    m:=BoundaryMult(n,AbsInt(F[i][1]),F[i][2]);
  b[i]:=Set(List([1..Length(m)],k->[AbsInt(m[k][1]),m[k][2]]));
 od;
 count:=0;
 lst:=Union(lst,b[1]);
 flag[1]:=1;
 count:=1;
 d:=1;
 while count<Length(F) do
  for cnr in [2..Length(F)] do
   if flag[cnr]=0 then
    if not IsEmpty(Intersection(lst,b[cnr])) then
     flag[cnr]:=1;
     count:=count+1;
     lst:=Union(lst,b[cnr]);
     break;
    fi;
   fi;
  od;
  d:=d+1;

  if not count=d then return false; fi;
 od;
 return true;
end;


############################ END CHECK##########################

    ############################################################################
    # Sort the (m-1)-faces of sigma into orbits under the action of Gamma_sigma


        Orb:=[];
        Orb[1]:=[];
        OrbElm:=[];
        OrbElm[1]:=[];
 ## Start with the cells of boundary(sigma) = bdry[1] :
        Add(Orb[1],bdry[1]);
        Add(OrbElm[1],id);

        for j in [2..Length(bdry)] do
            t:=0;
            for j1 in [1..Length(Orb)] do
                s:=IsSameOrbit(bdry[j],Orb[j1][1]);
                if not (s=false) then
                    Add(Orb[j1],bdry[j]);
                    Add(OrbElm[j1],s);
                    break;
                fi;
                t:=t+1;
            od;
            if t=Length(Orb) then Add(Orb,[bdry[j]]);Add(OrbElm,[id]);fi;
        od;


     ##################################################################
     # "Divide the boundary into rigid parts in the big cell [k,i]"
     #
     # Choose a fundamental domain T for the boundary of sigma
     # under the action of Gamma_sigma,
     # and tessellate the boundary of sigma with copies of T.
     ##################################################################

        Flag:=[];
        for j in [1..Length(Orb)] do
            Flag[j]:=[];
            for j1 in [1..Length(Orb[j])] do
                Flag[j][j1]:=0;
            od;
        od;

        for j in [1..Length(Orb[1])] do
            Flag[1][j]:=1;
        od;

        OrbFlag:=[];
        for j1 in [1..Length(Orb)] do
               OrbFlag[j1]:=0;
        od;
        OrbFlag[1]:=1;
if (k<=N) then
        T :=[Orb[1][1]];
        # In each orbit {g_{jt} sigma_j with j >= 2,
 # choose one cell such that its union with the cells in T is connected,
 # and such that the naive Euler characteristic of this union remains 1.



        T := SearchComponent(T[1]);


        if not Length(T)=Length(Orb) then
            Print("    Can not find a contractible fundamental domain! \n");
            Print(1/0);

        fi;


     ###################################################################
        w:=[];
        BoundaryTemp:=[];
        DimTemp:=[];
        for j in [0..(k)] do
            BoundaryTemp[j+1]:=[];
            DimTemp[j+1]:=0;
        od;

 # Use Algorithm 3 to construct a fundamental domain F for sigma


        # Create a list of already-taken-envelop
        Tcells:=[];
        for j in [1..k] do
          Tcells[j]:=[];
        od;
        convexTcells:=[];
        for j in [1..k] do
          convexTcells[j]:=[];
        od;

        d:=ConnectToCenter([k-1,T]);

        for j in [1..Length(OrbElm[1])] do
            Add(w,[SignInt(OrbElm[1][j])*d[1],pos(Elts[AbsInt(OrbElm[1][j])]*Elts[d[2]])]);
        od;
        return w;
      else

        T :=List(Orb,i->i[1]);
     ###################################################################
        w:=[];
        BoundaryTemp:=[];
        DimTemp:=[];
        for j in [0..(k)] do
            BoundaryTemp[j+1]:=[];
            DimTemp[j+1]:=0;
        od;

  # Use Algorithm 3 to construct a fundamental domain F for sigma
        F:=[];
        Tcells:=[];
        for j in [1..k] do
          Tcells[j]:=[];
        od;
        convexTcells:=[];
        for j in [1..k] do
          convexTcells[j]:=[];
        od;
        for tau in T do
             Ctau:=ConnectToCenter([k-1,[tau]]);
             Add(F,Ctau);
        od;

        for j in [1..Length(OrbElm[1])] do
            Append(w,List(F,d->[SignInt(OrbElm[1][j])*d[1],pos(Elts[AbsInt(OrbElm[1][j])]*Elts[d[2]])]));
        od;

        return w;
      fi;
    end;
    ##################################################################
    # Replacing a cell by its subdivision
    ReplaceCell:=function(k,m)
    local i, j, p, w, x, bdry, y, ww;
        w:=ShallowCopy(SubdividingCell(k,m));
        if k=N then
            Partition:=StructuralCopy(w);
            for i in [1..Length(Partition)] do
                Partition[i][1]:=AbsInt(Partition[i][1]-SignInt(Partition[i][1]));
            od;
        else Partition:=fail;
        fi;

        if k<=M and k<N then
        for i in [1..DimRec[k+2]] do
            bdry:=ShallowCopy(BoundaryRec[k+1][i]);
            p:=PositionsProperty(bdry,w->AbsInt(w[1])=m);
            for j in p do
                x:=bdry[j];
                ww:=ShallowCopy(w);
                if x[1]<0 then ww:=NegateWord(ww);fi;
                ww:=Mult(ww,k,x[2]);
                Append(bdry,ww);
            od;
            y:=bdry{p};
            bdry:=Set(bdry);
            SubtractSet(bdry,y);
            BoundaryRec[k+1][i]:=bdry;
        od;
        fi;
        BoundaryRec[k][m]:="del";
        StabRec[k+1][m]:="del";
    end;
    ##################################################################
    # Calculate the naive Euler characteristic for the list of k-cells L
    NaiveEulerChar:=function(k,L)
    local Cells, nrCells,x,w,j,s,id, Mat,t,M,b,p,h,r;

    Cells:=[];
    # id:=pos(One(C!.group));
    for j in [1..k+1] do
        Cells[j]:=[];

    od;

    j:=k;
    Append(Cells[k+1],L);


    while j>0 do
        for s in [1..Length(Cells[j+1])] do
            x:=Cells[j+1][s];
            w:=StructuralCopy(BoundaryRec[j][AbsInt(x[1])]);
            w:=Mult(w,j-1,x[2]);
            w:=List(w,a->[AbsInt(a[1]),a[2]]);
            Cells[j]:=Union(Cells[j],w);
        od;
        j:=j-1;
    od;
    nrCells:=List([1..k+1],a->Length(Cells[a]));
    w:=0;
    for j in [1..Length(nrCells)] do
        w:=w+(-1)^(j-1)*nrCells[j];
    od;

     return [w,nrCells];


    end;
    ################
    IsContractible:=function(k,L)
      local Boundaries,Coboundaries, Cells, x, w, b, j, s, t, Y, crit, l, flag,
            count, d, S, nrCells, eulerChar;



      # Construct the list of cells and the corresponding boundary and
      # coboundary of those cells
      Cells:=[];

      for j in [1..k+1] do
          Cells[j]:=[];

      od;
      j:=k;
      Append(Cells[k+1],L);
      Boundaries:=[];

      while j>0 do
          w:=[];
          for s in [1..Length(Cells[j+1])] do
              x:=Cells[j+1][s];
              w[s]:=BoundaryMult(j,AbsInt(x[1]),x[2]);
              w[s]:=List(w[s],a->[AbsInt(a[1]),a[2]]);
              Cells[j]:=Union(Cells[j],w[s]);
          od;
          if j=k then
              flag:=[];
              l:=Length(Cells[j+1]);
              for s in [1..l] do
                flag[s]:=0;
              od;
              flag[1]:=1;
              count:=1;
              S:=w[1];
              d:=1;
              while count<l do
                for s in [2..l] do
                  if flag[s]=0 then
                    if not IsEmpty(Intersection(S,w[s])) then
                      flag[s]:=1;
                      count:=count+1;
                      S:=Union(S,w[s]);
                      break;
                    fi;
                  fi;
                od;
                d:=d+1;
                if not count=d then return false;fi;
              od;
          fi;

          Boundaries[j+1]:=[];
          for s in [1..Length(Cells[j+1])] do
                  b:=List(w[s],a->Position(Cells[j],a));
                  Boundaries[j+1][s]:=Concatenation([Length(b)],b);
          od;

          j:=j-1;
      od;
      nrCells:=List([1..k+1],i->Length(Cells[i]));
      # Naive Euler Char of L=Union of T and \tau
      eulerChar:=Sum(List([1..k+1],i->(-1)^(i+1)*nrCells[i]));
      if not eulerChar=1 then return false;fi;

      Boundaries[1]:=List(Cells[1],a->[1,0]); # denote the boundary of 0-cell
      Boundaries[k+2]:=[];
      ######BOUNDARIES END########
      ######COBOUNDARIES BEGIN####
      Coboundaries:=[];; #Coboundaries[n+1] contains the info on n-cells.
      for j in [0..k] do
#        t:=j+3;
        Coboundaries[j+1]:=List(Boundaries[j+1],i->[0]);
        for s in [1..Length(Boundaries[j+2])] do
          b:=Boundaries[j+2][s];
          t:=Length(b);
          for i in b{[2..t]} do
            Coboundaries[j+1][i][1]:=Coboundaries[j+1][i][1]+1;
            Add(Coboundaries[j+1][i],s);
          od;
        od;

      od;
      Coboundaries[k+1]:=List(Boundaries[k+1],a->[0]);
      #####COBOUNDARIES END######
      Y:=ListsOfCellsToRegularCWComplex(k,Boundaries,Coboundaries,fail);
      crit:=CriticalCellsOfRegularCWComplex(Y);
      if (Length(crit)=1 and crit[1][1]=0) then return true;
      else return false;
      fi;

    end;
    ################
    IsSimplyConnected:=function(k,L,bool)
      local Boundaries,Coboundaries, Cells, x, w, b, j, s, t, Y, crit, l, flag,
            count, d, S, nrCells, eulerChar, P, Orientation;
      # Construct the list of cells and the corresponding boundary and
      # coboundary of those cells
      Cells:=[];

      for j in [1..k+1] do
          Cells[j]:=[];

      od;
      j:=k;
      Append(Cells[k+1],L);
      Boundaries:=[];
      Orientation:=[];

      while j>0 do
          w:=[];
          Orientation[j+1]:=[];
          for s in [1..Length(Cells[j+1])] do
              x:=Cells[j+1][s];
              w[s]:=BoundaryMult(j,AbsInt(x[1]),x[2]);
              Orientation[j+1][s]:=List(w[s],a->SignInt(a[1]));
              w[s]:=List(w[s],a->[AbsInt(a[1]),a[2]]);
              Cells[j]:=Union(Cells[j],w[s]);
          od;

          Boundaries[j+1]:=[];

          for s in [1..Length(Cells[j+1])] do
                  b:=List(w[s],a->Position(Cells[j],a));

                  Boundaries[j+1][s]:=Concatenation([Length(b)],b);
          od;


          j:=j-1;
      od;
      nrCells:=List([1..k+1],i->Length(Cells[i]));
      # Naive Euler Char of the boundary S of L
      eulerChar:=Sum(List([1..k+1],i->(-1)^(i+1)*nrCells[i]));
      if not eulerChar=1+(-1)^k then return false;fi;


      if bool=true then # bool=true if dimension of cells in S >=2 and
                        # Length(T)=Length(Orb)-1 so that the function is
                        # seeking the last cell in T


      Boundaries[1]:=List(Cells[1],a->[1,0]); # denote the boundary of 0-cell
      Orientation[1]:=List(Cells[1],a->[1]);
      Boundaries[k+2]:=[];
      ######BOUNDARIES END########
      ######COBOUNDARIES BEGIN####
      Coboundaries:=[];; #Coboundaries[n+1] contains the info on n-cells.
      for j in [0..k] do

        Coboundaries[j+1]:=List(Boundaries[j+1],i->[0]);
        for s in [1..Length(Boundaries[j+2])] do
          b:=Boundaries[j+2][s];
          t:=Length(b);
          for i in b{[2..t]} do
            Coboundaries[j+1][i][1]:=Coboundaries[j+1][i][1]+1;
            Add(Coboundaries[j+1][i],s);
          od;
        od;

      od;
      Coboundaries[k+1]:=List(Boundaries[k+1],a->[0]);
      #####COBOUNDARIES END######
      Y:=ListsOfCellsToRegularCWComplex(k,Boundaries,Coboundaries,Orientation);

      P:=FundamentalGroupOfRegularCWComplex(Y);

      if (IsTrivial(P)) then return true;
      else return false;
      fi;
    fi;
    return true;

    end;
    ###############
    CellHomology:=function(k,L)
    local Cells, nrCells,x,w,j,s,id, Mat,t,M,b,p,h,r;

    Cells:=[];

    for j in [1..k+1] do
        Cells[j]:=[];

    od;

    j:=k;
    Append(Cells[k+1],L);


    while j>0 do
        for s in [1..Length(Cells[j+1])] do
            x:=Cells[j+1][s];

            w[s]:=BoundaryMult(j,AbsInt(x[1]),x[2]);
            w:=List(w,a->[AbsInt(a[1]),a[2]]);
            Cells[j]:=Union(Cells[j],w);
        od;
        j:=j-1;
    od;
    nrCells:=List([1..k+1],a->Length(Cells[a]));

  Mat:=[];
  for j in [2..k+1] do
    M:=[];
    for s in [1..Length(Cells[j])] do
       M[s]:=[];
       for t in [1..Length(Cells[j-1])] do
          M[s][t]:=0;
       od;
    od;

    for s in [1..Length(Cells[j])] do
        x:=Cells[j][s];

        b:=BoundaryMult(j-1,AbsInt(x[1]),x[2]);
        for t in [1..Length(b)] do
            p:=Position(Cells[j-1],[AbsInt(b[t][1]),b[t][2]]);
            M[s][p]:=SignInt(b[t][1]);
        od;
    od;
    Mat[j]:=SmithNormalFormIntegerMat(TransposedMat(M));


  od;
    w:=[];
    # Create Mat[1] for d_0
    M:=[];
    for s in [1..Length(Cells[1])] do
        M[s]:=0;
    od;
    Mat[1]:=TransposedMat([M]);

    # Create Mat[k+1] do d_(k+1)
    M:=[];
    for s in [1..Length(Cells[k+1])] do
        M[s]:=0;
    od;
    Mat[k+2]:=[M];

    for j in [1..k+1] do
        h:=[];
        r:=Length(Cells[j])-RankMat(Mat[j])-RankMat(Mat[j+1]);
        for s in [1..r] do
     h[s]:=0;
        od;
        for s in [1..RankMat(Mat[j+1])] do
            if not Mat[j+1][s][s]=1 then
                Add(h,Mat[j+1][s][s]);
            fi;
        od;
        Add(w,h);

    od;



     return w;


    end;



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

    ##################################################################
    # Main part: subdividing the fundamental domain
    NotRigid:=[];
    dims:=ShallowCopy(DimRec);
    i:=1;

    while i<=M do

        j:=1;
        while j<=dims[i+1] do
            if not IsRigidCell(i,j) then
                Add(NotRigid,[i,j]);
            fi;
            j:=j+1;
        od;

        i:=i+1;

    od;
    for x in NotRigid do
#       Apply Rigid Facets Subdivision to the cell x :
        ReplaceCell(x[1],x[2]);
    od;


    #Delete cells which are already replaced by its subdivision

    t:=1;
    for w in [1..Length(NotRigid)] do
        k:=NotRigid[w][1];
        j:=NotRigid[w][2];
        if k<N then
        for i in [1..DimRec[k+2]] do
            bdry:=BoundaryRec[k+1][i];

            if not IsString(bdry) then
               for x in bdry do
                    if AbsInt(x[1])>j then
                        x[1]:=x[1]-SignInt(x[1]);
                    fi;
                od;
             fi;
             BoundaryRec[k+1][i]:=bdry;

         od;
         fi;
         dims[k+1]:=dims[k+1]-1;
         DimRec[k+1]:=DimRec[k+1]-1;
         Remove(BoundaryRec[k],j);
         Remove(StabRec[k+1],j);
         if IsBound(NotRigid[w+1]) and NotRigid[w+1][1]=NotRigid[w][1] then
             NotRigid[w+1][2]:=NotRigid[w+1][2]-t;
             t:=t+1;
         else
             t:=1;
         fi;

    od;

    ##################################################################
    Boundary:=function(k,m)
        if m<0 then return NegateWord(BoundaryRec[k][AbsInt(m)]);fi;
        return BoundaryRec[k][m];
    end;

    Stabilizer:=function(k,m)
        return StabRec[k+1][m];
    end;

    Dimension:=function(k)
        if k>N then return 0;fi;
        return DimRec[k+1];
    end;

    Action:=function(k,i,j)
        return 1;
    end;
    ##################################################################
return Objectify(HapNonFreeResolution,
    rec(
    dimension:=Dimension,
    Partition:=Partition,
    boundary:=Boundary,
    homotopy:=fail,
    elts:=Elts,
    group:=C!.group,
    stabilizer:=Stabilizer,
    action:=Action,
    subdividing:=SubdividingCell,
    replacecell:=ReplaceCell,
    isrigid:=IsRigidCell,
    properties:=
    [["length",Maximum(1000,N)],
    ["characteristic",0],
    ["type","resolution"]]  ));
end);


################### end of ControlledSubdivision ############################


################## HybridSubdivision #######################################
## Hybrid Subdivision combine the goods of both RFS and VSS. It will Use
## RFS to subdivide the cells in dimension N-1
## and continue using VSS to subdivide top dimensional cells.

DeclareGlobalFunction("HybridSubdivision");
InstallGlobalFunction(HybridSubdivision,
function(arg)
local W, StabRec, i, j, N, M, x, bdry, s1, s2, p, k, w, t,
      DimRec, BoundaryRec, id, dims, NotRigid, NewCell,
      Cell, Elts, Boundary, Dimension, CLeftCosetElt, LCoset,
      pos, Stab, Mult,DimTemp, BoundaryTemp, Partition, ConnectedCheck,
      Stabilizer, Action, IsRigidCell, ReplaceCell, SubdividingCell,
      Orbit, C, NaiveEulerChar, CellHomology, exportFile, IsContractible,
      CosetRec, IsSimplyConnected, BoundaryMult, BoundaryMultRec, jj, jjj;

    C:=arg[1];
    i:=0;
    while C!.dimension(i)>0 do
        i:=i+1;
    od;
    N:=i-1; # Length of the chain complex
    M:=N;

if Length(arg)=2 then M:=Minimum(N,arg[2]);
fi;


    Elts:=C!.elts;
    StabRec:=[];
    DimRec:=[];
    CosetRec:=[];

    ##################################################################
    # If g in Elts return the position of g in the list,
    # otherwise, add g to Elts and return the position.
    pos:=function(g)
    local posit;

        posit:=Position(Elts,g);
        if posit=fail then
            Add(Elts,g);
            return Length(Elts);
        else
            return posit;
        fi;
    end;
    ##################################################################
    id:=pos(One(C!.group));
    ##################################################################
    # return the stabilizer of g*e,
    #
    Stab:=function(e,g)
    return ConjugateGroup(StabRec[e[1]+1][AbsInt(e[2])],Elts[g]^-1);
    end;
    ##################################################################
    # returns  a  "canonical"  representative  of  the  right  coset
    # Elts[g]*Stab[i+1][j]

    ######Add a record for left coset
    for i in [0..N] do
      CosetRec[i+1]:=[];
    od;

    CLeftCosetElt:=function(i,j,g)
      local p;
      if IsBound(CosetRec[AbsInt(i)+1][j]) then
        if IsBound(CosetRec[AbsInt(i)+1][j][g]) then
          return CosetRec[AbsInt(i)+1][j][g];
        else
          p:=pos(CanonicalRightCountableCosetElement
                                  (StabRec[AbsInt(i)+1][j],Elts[g]^-1)^-1);
          CosetRec[AbsInt(i)+1][j][g]:=p;
          return p;
        fi;
      else
        CosetRec[i+1][j]:=[];
        p:=pos(CanonicalRightCountableCosetElement
                                (StabRec[AbsInt(i)+1][j],Elts[g]^-1)^-1);
        CosetRec[AbsInt(i)+1][j][g]:=p;
        return p;
      fi;

    end;
    ##################################################################
    ##
    ##  Input:  A list L, degree k, position g of an element
    ##  Output: Product of g and L.
    ##
    Mult:=function(L,k,g)
    local x,w,t,h,y,vv,LL;
        vv:=[];
        LL:=ShallowCopy(L);
        for x in [1..Length(LL)] do
            w:=Elts[g]*Elts[LL[x][2]];
            t:=CLeftCosetElt(k,AbsInt(LL[x][1]),pos(w));
            Add(vv,[LL[x][1],t]);
        od;
        return vv;
    end;
    ###################################################################
    # Store essential data: stabilizers, boundaries, dimensions

    Orbit:=[];
    for i in [1..N+1] do
        Orbit[i]:=[];
    od;

    for i in [0..N] do
        StabRec[i+1]:=[];
        DimRec[i+1]:=C!.dimension(i);
        for j in [1..C!.dimension(i)] do
            StabRec[i+1][j]:=C!.stabilizer(i,j);
        od;
    od;

    BoundaryRec:=[];
    for i in [1..N] do
        BoundaryRec[i]:=[];
        for j in [1..DimRec[i+1]] do
            bdry:=C!.boundary(i,j);
            BoundaryRec[i][j]:=[];
            for x in bdry do
                s1:=C!.action(i-1,AbsInt(x[1]),x[2]);
                p:=pos(CanonicalRightCountableCosetElement
                            (C!.stabilizer(i-1,AbsInt(x[1])),Elts[x[2]]^-1)^-1);
                s2:=C!.action(i-1,AbsInt(x[1]),p);

                Add(BoundaryRec[i][j],[s1*s2*x[1],p]);
            od;

        od;
    od;

    #############
    BoundaryMultRec:=[];
    for jj in [1..N] do
      BoundaryMultRec[jj]:=[];
      for jjj in [1..Length(BoundaryRec[jj])] do
        BoundaryMultRec[jj][jjj]:=[];
      od;
    od;
    #############
    BoundaryMult:=function(i,j,g)
      local w;
      if IsBound(BoundaryMultRec[i]) then
        if IsBound(BoundaryMultRec[i][j]) then
          if IsBound(BoundaryMultRec[i][j][g]) then
            return BoundaryMultRec[i][j][g];
          else
            BoundaryMultRec[i][j][g]:=Mult(BoundaryRec[i][j],i-1,g);
            return BoundaryMultRec[i][j][g];
          fi;
        else
          BoundaryMultRec[i][j]:=[];
          BoundaryMultRec[i][j][g]:=Mult(BoundaryRec[i][j],i-1,g);
          return BoundaryMultRec[i][j][g];
        fi;
      else
        BoundaryMultRec[i]:=[];
        BoundaryMultRec[i][j]:=[];
        BoundaryMultRec[i][j][g]:=Mult(BoundaryRec[i][j],g);
        return BoundaryMultRec[i][j][g];
      fi;
    end;
    ##################################################################


    ##################################################################
    # Check whether the cell is rigid or not

    IsRigidCell:=function(k,m)
    local bdry, intst, L;
        bdry:=BoundaryRec[k][m];
        L:=List(bdry,w->Elements(ConjugateGroup(StabRec[k][AbsInt(w[1])],Elts[w[2]]^-1)));
        intst:=Intersection(L);
        if not Elements(StabRec[k+1][m])=Elements(intst) then
            return false;
        else return true;
        fi;

    end;
    ##################################################################


    ##################################################################
    # Subdividing (k+1)-cell number i (denoted by sigma) with the
    # Rigid Facets Subdivision Algorithm :
    SubdividingCell:=function(k,i)
    local bdry, w, x, d, y, a, b, T, j, j1,j2, z, t, c, Flag, s, s1, ConnectToCenter,
          SearchComponent,temp, l1, bdry1, bdry2, Orb, IsSameOrbit, OrbFlag, OrbElm,
          ConnectedCheck,tau,Ctau,F,exportBoundary,exportStabilizer,
          Boundaries,Stabilizers,Tcells, convexTcells,p;

 # Record the barycenter of sigma as a vertex, with stabilizer Gamma_sigma :
        DimRec[1]:=DimRec[1]+1;
        Add(StabRec[1],StabRec[k+1][i]);

        bdry:=ShallowCopy(BoundaryRec[k][i]);


        #############################################################
        # In each orbit {g_{jt} sigma_j with j >= 2,
 # choose one cell such that its union with the cells in T is connected,
 # and such that the naive Euler characteristic of this union remains 1.
  SearchComponent:=function(q)
  local x,z,j,w,TT, tau, boundaryOFtau, S, Ps, ps,countPrint, bool;
      x:=[];
      w:=[];
      Ps:=[[1,1]];

      j:=1;
      while not Length(T)=Length(Orb) do


          while j<=Length(Orb) do
              if OrbFlag[j]=0 then

                  while j1<=Length(Orb[j]) do
                      if Flag[j][j1]=0 then
                          TT:=Union(T,[Orb[j][j1]]);



                          if IsContractible(k-1,TT) then

                            if k>=2 then
                              S:=[];
                              for tau in TT do

                                   boundaryOFtau:=BoundaryMult(k-1,AbsInt(tau[1]),tau[2]);
                                if  tau[1]<0 then
                                    boundaryOFtau:=NegateWord(boundaryOFtau);
                                fi;
                                Append(S,boundaryOFtau);

                              od;
                              S:=AlgebraicReduction(S);
                              if (Length(T)=Length(Orb)-1) and k>=4 then bool:=true;
                              else bool:=false;
                              fi;
                              if IsSimplyConnected(k-2,S,bool) then

                                  Add(T,Orb[j][j1]);
                                  Flag[j][j1]:=1;
                                  Add(Ps,[j,j1]);
                                  OrbFlag[j]:=1;
                                  j:=1;
                                  break;
                              fi;
                            else
                              Add(T,Orb[j][j1]);
                              Flag[j][j1]:=1;
                              Add(Ps,[j,j1]);
                              OrbFlag[j]:=1;
                              j:=1;
                              break;

                            fi;

                          fi;

                      fi;

                      j1:=j1+1;
                    od;
                  if j1<=Length(Orb[j]) then
                      break;

                  fi;
              fi;
              j:=j+1;
              j1:=1;
          od;
          if j>Length(Orb) then
              Remove(T);

              ps:=Ps[Length(Ps)];
              OrbFlag[ps[1]]:=0;
              Flag[ps[1]][ps[2]]:=0;
              Remove(Ps);
              j:=ps[1];
              j1:=ps[2]+1;
              while j1>Length(Orb[j]) do
                ps:=Ps[Length(Ps)];
                OrbFlag[ps[1]]:=0;
                Flag[ps[1]][ps[2]]:=0;
                Remove(Ps);
                Remove(T);
                j:=ps[1];
                j1:=ps[2]+1;

              od;
              if Length(T)=0 then
                  Print("\n Could not find any fundamental domain for the ",i,"th of the",k,"-cells! \n");
                  Print(1/0);
              fi;
          fi;

      od;

      return T;
  end;
        #############################################################

    #################################################################
    # Connect the collection of (k-1)-cells T to the barycenter of the cell sigma.
    # The elements of T, as well as sigma, are in the format [k,i,gamma]:
    # dimension k, obtained by sending the ith orbit representative to its image under the element gamma (specified by its number in Gamma via C!.elts).
    # Input: e := [k-1, T].
    ConnectToCenter:=function(e)
    local bdry1, x, stab, bdrye, w, stablst, redbdry, tau, boundaryOFtau,
          LCoset, AddCell;


    ##################################################################
    # Add a k-cell with stabilizer stab and boundary bdry
    # to the cell complex
    AddCell:=function(m,stab,bdry,e)
    local a,j,g,s,w,u,v;

    if m<k then
        w:=e[1];

        for j in [1..DimTemp[m+1]] do
            u:=BoundaryTemp[m+1][j];
     # The cells u and w are in the format [orbit number, gamma],
     # such that gamma (pointing to C!.elts) sends the orbit representative to the cell.
     # We make use of the common vertex, namely the barycenter of sigma,
     # to test if there exists gamma in Gamma_sigma with gamma u = w.
            s:=DimRec[m+1]-DimTemp[m+1]+j;
            if AbsInt(w[1])=AbsInt(u[1]) then


              for v in StabRec[m][AbsInt(u[1])] do
                g:=Elts[w[2]]*v*Elts[u[2]]^-1;
                 if g in StabRec[k+1][i] then
                    return [s, CLeftCosetElt(m,s,pos(g))];
                  fi;
              od;

            fi;
        od;

        DimTemp[m+1]:=DimTemp[m+1]+1;

        Add(BoundaryTemp[m+1],e[1]);

 #Update the dimension :
        DimRec[m+1]:=DimRec[m+1]+1;

 # Record the stabilizer :
        Add(StabRec[m+1],stab);\

 # Record the boundary of F :
        Add(BoundaryRec[m],bdry);

    else
        DimRec[m+1]:=DimRec[m+1]+1;
        Add(StabRec[m+1],stab);
         Add(BoundaryRec[m],bdry);
    fi;
    return [DimRec[m+1],CLeftCosetElt(m,DimRec[m+1],id)];
    end;
    ##################################################################

        if e[1]=0 # k-1 = 0 : e[2] = T is a collection of vertices.
 then
            p:=Position(Tcells[1],e[2][1]);
            if not p=fail then return convexTcells[1][p];
            else
            bdry1:=[[-SignInt(e[2][1][1])*DimRec[1],CLeftCosetElt(0,DimRec[1],id)],[e[2][1][1],e[2][1][2]]];

            stab:=Intersection(Stab([0,e[2][1][1]],e[2][1][2]),StabRec[k+1][i]);
            w:=AddCell(1,stab,bdry1,e[2]);
            Add(Tcells[1],e[2][1]);
            Add(convexTcells[1],w);
            return w;
            fi;
        fi;

        stablst:=[];
        bdry1:=[];
        Append(bdry1,e[2]);
        bdrye:=[];

 ## Enumerate the set S := boundary(geometric realization of T) =: bdrye,
 ## by letting tau run through the cells of T = e[2],
 ## and afterwards cancelling cells in S which appear in pairs with opposite orientation.
        for tau in e[2] do
            # boundaryOFtau:=Mult(BoundaryRec[e[1]][AbsInt(tau[1])],e[1]-1,tau[2]);
            boundaryOFtau:=BoundaryMult(e[1],AbsInt(tau[1]),tau[2]);
            if  tau[1]<0 then boundaryOFtau:=NegateWord(boundaryOFtau);fi;
            Append(bdrye,boundaryOFtau);
            Add(stablst,Stab([e[1],tau[1]],tau[2]));
        od;
        Add(stablst,StabRec[k+1][i]);
 ## Cancel cells in S which appear in pairs with opposite orientation :
        bdrye:=AlgebraicReduction(bdrye);
 ## Construction of S = bdrye is now completed.

 ## For every x in S, take the convex envelope w := e(x) of x and the barycenter of sigma :

        for x in bdrye do
          p:=Position(Tcells[e[1]],[AbsInt(x[1]),x[2]]);
          if not p=fail then
            w:=convexTcells[e[1]][p];
          else
            w:=ConnectToCenter([e[1]-1,[[AbsInt(x[1]),x[2]]]]);
            Add(Tcells[e[1]],[AbsInt(x[1]),x[2]]);
            Add(convexTcells[e[1]],w);
          fi;
            Add(bdry1,[-SignInt(x[1])*w[1],w[2]]);
     ## Here, the orientation of w := e(x) has been recorded in -SignInt(x[1]).
        od;


 # Calculate the stabilizer :
        stab:=Intersection(stablst);
        if not Length(e[2])=1 then

            w:=AddCell(e[1]+1,stab,bdry1,e[2]);

            return w;
        else
            return AddCell(e[1]+1,stab,bdry1,e[2]);
        fi;
    end;
    ##################################################################
    ## End of function ConnectToCenter, end of Algorithm 3.
    ##################################################################


    IsSameOrbit:=function(e,f)
    local s1,s;
        if AbsInt(e[1])=AbsInt(f[1]) then
            s:=StabRec[k+1][i];
            s1:=StabRec[k][AbsInt(e[1])];
            for x in s do
                if Elts[f[2]]^-1*x*Elts[e[2]] in s1 then
                    return SignInt(e[1])*SignInt(f[1])*pos(x);
                fi;
            od;
            return false;
         else
            return false;
         fi;
    end;

##### 2017 - LUXEMBOURG


##### Check if the fundamental domain formed by sub[1] is connected ####

ConnectedCheck:=function(F,n)
 local b, lst, count, flag, i, d, cnr, m;

 ###

 lst:=[];
        flag:=[];
 b:=[];
 for i in [1..Length(F)] do
  flag[i]:=0;

    m:=BoundaryMult(n,AbsInt(F[i][1]),F[i][2]);
  b[i]:=Set(List([1..Length(m)],k->[AbsInt(m[k][1]),m[k][2]]));
 od;
 count:=0;
 lst:=Union(lst,b[1]);
 flag[1]:=1;
 count:=1;
 d:=1;
 while count<Length(F) do
  for cnr in [2..Length(F)] do
   if flag[cnr]=0 then
    if not IsEmpty(Intersection(lst,b[cnr])) then
     flag[cnr]:=1;
     count:=count+1;
     lst:=Union(lst,b[cnr]);
     break;
    fi;
   fi;
  od;
  d:=d+1;

  if not count=d then return false; fi;
 od;
 return true;
end;


############################ END CHECK##########################

    ############################################################################
    # Sort the (m-1)-faces of sigma into orbits under the action of Gamma_sigma
        Orb:=[];
        Orb[1]:=[];
        OrbElm:=[];
        OrbElm[1]:=[];
 ## Start with the cells of boundary(sigma) = bdry[1] :
        Add(Orb[1],bdry[1]);
        Add(OrbElm[1],id);

        for j in [2..Length(bdry)] do
            t:=0;
            for j1 in [1..Length(Orb)] do
                s:=IsSameOrbit(bdry[j],Orb[j1][1]);
                if not (s=false) then
                    Add(Orb[j1],bdry[j]);
                    Add(OrbElm[j1],s);
                    break;
                fi;
                t:=t+1;
            od;
            if t=Length(Orb) then Add(Orb,[bdry[j]]);Add(OrbElm,[id]);fi;
        od;

     ##################################################################
     # "Divide the boundary into rigid parts in the big cell [k,i]"
     #
     # Choose a fundamental domain T for the boundary of sigma
     # under the action of Gamma_sigma,
     # and tessellate the boundary of sigma with copies of T.
     ##################################################################

        Flag:=[];
        for j in [1..Length(Orb)] do
            Flag[j]:=[];
            for j1 in [1..Length(Orb[j])] do
                Flag[j][j1]:=0;
            od;
        od;

        for j in [1..Length(Orb[1])] do
            Flag[1][j]:=1;
        od;

        OrbFlag:=[];
        for j1 in [1..Length(Orb)] do
               OrbFlag[j1]:=0;
        od;
        OrbFlag[1]:=1;
if (k<=N-1) then

        T :=[Orb[1][1]];
        # In each orbit {g_{jt} sigma_j with j >= 2,
 # choose one cell such that its union with the cells in T is connected,
 # and such that the naive Euler characteristic of this union remains 1.

        T := SearchComponent(T[1]);


        if not Length(T)=Length(Orb) then
            Print("    Can not find a contractible fundamental domain! \n");
            Print(1/0);
        else


        fi;


     ###################################################################
        w:=[];
        BoundaryTemp:=[];
        DimTemp:=[];
        for j in [0..(k)] do
            BoundaryTemp[j+1]:=[];
            DimTemp[j+1]:=0;
        od;

 # Use Algorithm 3 to construct a fundamental domain F for sigma

        # Create a list of already-taken-envelop
        Tcells:=[];
        for j in [1..k] do
          Tcells[j]:=[];
        od;
        convexTcells:=[];
        for j in [1..k] do
          convexTcells[j]:=[];
        od;

        d:=ConnectToCenter([k-1,T]);

        for j in [1..Length(OrbElm[1])] do
            Add(w,[SignInt(OrbElm[1][j])*d[1],pos(Elts[AbsInt(OrbElm[1][j])]*Elts[d[2]])]);
        od;
        return w;
      else

        T :=List(Orb,i->i[1]);
     ###################################################################
        w:=[];
        BoundaryTemp:=[];
        DimTemp:=[];
        for j in [0..(k)] do
            BoundaryTemp[j+1]:=[];
            DimTemp[j+1]:=0;
        od;

  # Use Algorithm 3 to construct a fundamental domain F for sigma
        F:=[];
        Tcells:=[];
        for j in [1..k] do
          Tcells[j]:=[];
        od;
        convexTcells:=[];
        for j in [1..k] do
          convexTcells[j]:=[];
        od;
        for tau in T do
             Ctau:=ConnectToCenter([k-1,[tau]]);
             Add(F,Ctau);
        od;

        for j in [1..Length(OrbElm[1])] do
            Append(w,List(F,d->[SignInt(OrbElm[1][j])*d[1],pos(Elts[AbsInt(OrbElm[1][j])]*Elts[d[2]])]));
        od;

        return w;
      fi;
    end;
    ##################################################################
    # Replacing a cell by its subdivision
    ReplaceCell:=function(k,m)
    local i, j, p, w, x, bdry, y, ww;
        w:=ShallowCopy(SubdividingCell(k,m));
        if k=N then
            Partition:=StructuralCopy(w);
            for i in [1..Length(Partition)] do
                Partition[i][1]:=AbsInt(Partition[i][1]-SignInt(Partition[i][1]));
            od;
        else Partition:=fail;
        fi;

        if k<=M and k<N then
        for i in [1..DimRec[k+2]] do
            bdry:=ShallowCopy(BoundaryRec[k+1][i]);
            p:=PositionsProperty(bdry,w->AbsInt(w[1])=m);
            for j in p do
                x:=bdry[j];
                ww:=ShallowCopy(w);
                if x[1]<0 then ww:=NegateWord(ww);fi;
                ww:=Mult(ww,k,x[2]);
                Append(bdry,ww);
            od;
            y:=bdry{p};
            bdry:=Set(bdry);
            SubtractSet(bdry,y);
            BoundaryRec[k+1][i]:=bdry;
        od;
        fi;
        BoundaryRec[k][m]:="del";
        StabRec[k+1][m]:="del";
    end;
    ##################################################################
    # Calculate the naive Euler characteristic for the list of k-cells L
    NaiveEulerChar:=function(k,L)
    local Cells, nrCells,x,w,j,s,id, Mat,t,M,b,p,h,r;

    Cells:=[];

    for j in [1..k+1] do
        Cells[j]:=[];

    od;

    j:=k;
    Append(Cells[k+1],L);


    while j>0 do
        for s in [1..Length(Cells[j+1])] do
            x:=Cells[j+1][s];
            w:=StructuralCopy(BoundaryRec[j][AbsInt(x[1])]);
            w:=Mult(w,j-1,x[2]);
            w:=List(w,a->[AbsInt(a[1]),a[2]]);
            Cells[j]:=Union(Cells[j],w);
        od;
        j:=j-1;
    od;
    nrCells:=List([1..k+1],a->Length(Cells[a]));
    w:=0;
    for j in [1..Length(nrCells)] do
        w:=w+(-1)^(j-1)*nrCells[j];
    od;



     return [w,nrCells];


    end;
    ################
    IsContractible:=function(k,L)
      local Boundaries,Coboundaries, Cells, x, w, b, j, s, t, Y, crit, l, flag,
            count, d, S, nrCells, eulerChar;



      # Construct the list of cells and the corresponding boundary and
      # coboundary of those cells
      Cells:=[];

      for j in [1..k+1] do
          Cells[j]:=[];

      od;
      j:=k;
      Append(Cells[k+1],L);
      Boundaries:=[];

      while j>0 do
          w:=[];
          for s in [1..Length(Cells[j+1])] do
              x:=Cells[j+1][s];

              w[s]:=BoundaryMult(j,AbsInt(x[1]),x[2]);
              w[s]:=List(w[s],a->[AbsInt(a[1]),a[2]]);
              Cells[j]:=Union(Cells[j],w[s]);
          od;
          if j=k then
              flag:=[];
              l:=Length(Cells[j+1]);
              for s in [1..l] do
                flag[s]:=0;
              od;
              flag[1]:=1;
              count:=1;
              S:=w[1];
              d:=1;
              while count<l do
                for s in [2..l] do
                  if flag[s]=0 then
                    if not IsEmpty(Intersection(S,w[s])) then
                      flag[s]:=1;
                      count:=count+1;
                      S:=Union(S,w[s]);
                      break;
                    fi;
                  fi;
                od;
                d:=d+1;
                if not count=d then return false;fi;
              od;
          fi;

          Boundaries[j+1]:=[];
          for s in [1..Length(Cells[j+1])] do
                  b:=List(w[s],a->Position(Cells[j],a));
                  Boundaries[j+1][s]:=Concatenation([Length(b)],b);
          od;


          j:=j-1;
      od;
      nrCells:=List([1..k+1],i->Length(Cells[i]));
      # Naive Euler Char of L=Union of T and \tau
      eulerChar:=Sum(List([1..k+1],i->(-1)^(i+1)*nrCells[i]));
      if not eulerChar=1 then return false;fi;
      #
      # # Naive Euler Char of the boundary S of L

      Boundaries[1]:=List(Cells[1],a->[1,0]); # denote the boundary of 0-cell
      Boundaries[k+2]:=[];
      ######BOUNDARIES END########
      ######COBOUNDARIES BEGIN####
      Coboundaries:=[];; #Coboundaries[n+1] contains the info on n-cells.
      for j in [0..k] do

        Coboundaries[j+1]:=List(Boundaries[j+1],i->[0]);
        for s in [1..Length(Boundaries[j+2])] do
          b:=Boundaries[j+2][s];
          t:=Length(b);
          for i in b{[2..t]} do
            Coboundaries[j+1][i][1]:=Coboundaries[j+1][i][1]+1;
            Add(Coboundaries[j+1][i],s);
          od;
        od;

      od;
      Coboundaries[k+1]:=List(Boundaries[k+1],a->[0]);
      #####COBOUNDARIES END######
      Y:=ListsOfCellsToRegularCWComplex(k,Boundaries,Coboundaries,fail);
      crit:=CriticalCellsOfRegularCWComplex(Y);

      if (Length(crit)=1 and crit[1][1]=0) then return true;
      else return false;
      fi;

    end;
    ################
    IsSimplyConnected:=function(k,L,bool)
      local Boundaries,Coboundaries, Cells, x, w, b, j, s, t, Y, crit, l, flag,
            count, d, S, nrCells, eulerChar, P, Orientation;



      # Construct the list of cells and the corresponding boundary and
      # coboundary of those cells
      Cells:=[];

      for j in [1..k+1] do
          Cells[j]:=[];

      od;
      j:=k;
      Append(Cells[k+1],L);
      Boundaries:=[];
      Orientation:=[];

      while j>0 do
          w:=[];
          Orientation[j+1]:=[];
          for s in [1..Length(Cells[j+1])] do
              x:=Cells[j+1][s];

              w[s]:=BoundaryMult(j,AbsInt(x[1]),x[2]);
              Orientation[j+1][s]:=List(w[s],a->SignInt(a[1]));
              w[s]:=List(w[s],a->[AbsInt(a[1]),a[2]]);
              Cells[j]:=Union(Cells[j],w[s]);
          od;

          Boundaries[j+1]:=[];

          for s in [1..Length(Cells[j+1])] do
                  b:=List(w[s],a->Position(Cells[j],a));

                  Boundaries[j+1][s]:=Concatenation([Length(b)],b);
          od;


          j:=j-1;
      od;
      nrCells:=List([1..k+1],i->Length(Cells[i]));

      eulerChar:=Sum(List([1..k+1],i->(-1)^(i+1)*nrCells[i]));
      if not eulerChar=1+(-1)^k then return false;fi;


      if bool=true then


      Boundaries[1]:=List(Cells[1],a->[1,0]); # denote the boundary of 0-cell
      Orientation[1]:=List(Cells[1],a->[1]);
      Boundaries[k+2]:=[];
      ######BOUNDARIES END########
      ######COBOUNDARIES BEGIN####
      Coboundaries:=[];; #Coboundaries[n+1] contains the info on n-cells.
      for j in [0..k] do

        Coboundaries[j+1]:=List(Boundaries[j+1],i->[0]);
        for s in [1..Length(Boundaries[j+2])] do
          b:=Boundaries[j+2][s];
          t:=Length(b);
          for i in b{[2..t]} do
            Coboundaries[j+1][i][1]:=Coboundaries[j+1][i][1]+1;
            Add(Coboundaries[j+1][i],s);
          od;
        od;

      od;
      Coboundaries[k+1]:=List(Boundaries[k+1],a->[0]);
      #####COBOUNDARIES END######
      Y:=ListsOfCellsToRegularCWComplex(k,Boundaries,Coboundaries,Orientation);

      P:=FundamentalGroupOfRegularCWComplex(Y);

      if (IsTrivial(P)) then return true;
      else return false;
      fi;
      fi;
      return true;
    end;
    ###############
    CellHomology:=function(k,L)
    local Cells, nrCells,x,w,j,s,id, Mat,t,M,b,p,h,r;

    Cells:=[];

    for j in [1..k+1] do
        Cells[j]:=[];

    od;

    j:=k;
    Append(Cells[k+1],L);


    while j>0 do
        for s in [1..Length(Cells[j+1])] do
            x:=Cells[j+1][s];

            w[s]:=BoundaryMult(j,AbsInt(x[1]),x[2]);
            w:=List(w,a->[AbsInt(a[1]),a[2]]);
            Cells[j]:=Union(Cells[j],w);
        od;
        j:=j-1;
    od;
    nrCells:=List([1..k+1],a->Length(Cells[a]));

  Mat:=[];
  for j in [2..k+1] do
    M:=[];
    for s in [1..Length(Cells[j])] do
       M[s]:=[];
       for t in [1..Length(Cells[j-1])] do
          M[s][t]:=0;
       od;
    od;

    for s in [1..Length(Cells[j])] do
        x:=Cells[j][s];

        b:=BoundaryMult(j-1,AbsInt(x[1]),x[2]);
        for t in [1..Length(b)] do
            p:=Position(Cells[j-1],[AbsInt(b[t][1]),b[t][2]]);
            M[s][p]:=SignInt(b[t][1]);
        od;
    od;
    Mat[j]:=SmithNormalFormIntegerMat(TransposedMat(M));


  od;
    w:=[];
    # Create Mat[1] for d_0
    M:=[];
    for s in [1..Length(Cells[1])] do
        M[s]:=0;
    od;
    Mat[1]:=TransposedMat([M]);

    # Create Mat[k+1] do d_(k+1)
    M:=[];
    for s in [1..Length(Cells[k+1])] do
        M[s]:=0;
    od;
    Mat[k+2]:=[M];

    for j in [1..k+1] do
        h:=[];
        r:=Length(Cells[j])-RankMat(Mat[j])-RankMat(Mat[j+1]);
        for s in [1..r] do
     h[s]:=0;
        od;
        for s in [1..RankMat(Mat[j+1])] do
            if not Mat[j+1][s][s]=1 then
                Add(h,Mat[j+1][s][s]);
            fi;
        od;
        Add(w,h);

    od;



     return w;


    end;



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

    ##################################################################
    # Main part: subdividing the fundamental domain
    NotRigid:=[];
    dims:=ShallowCopy(DimRec);
    i:=1;

    while i<=M do

        j:=1;
        while j<=dims[i+1] do
            if not IsRigidCell(i,j) then
                Add(NotRigid,[i,j]);
            fi;
            j:=j+1;
        od;

        i:=i+1;

    od;
    for x in NotRigid do
#       Apply Rigid Facets Subdivision to the cell x :
        ReplaceCell(x[1],x[2]);
    od;


    #Delete cells which are already replaced by its subdivision

    t:=1;
    for w in [1..Length(NotRigid)] do
        k:=NotRigid[w][1];
        j:=NotRigid[w][2];
        if k<N then
        for i in [1..DimRec[k+2]] do
            bdry:=BoundaryRec[k+1][i];

            if not IsString(bdry) then
               for x in bdry do
                    if AbsInt(x[1])>j then
                        x[1]:=x[1]-SignInt(x[1]);
                    fi;
                od;
             fi;
             BoundaryRec[k+1][i]:=bdry;

         od;
         fi;
         dims[k+1]:=dims[k+1]-1;
         DimRec[k+1]:=DimRec[k+1]-1;
         Remove(BoundaryRec[k],j);
         Remove(StabRec[k+1],j);
         if IsBound(NotRigid[w+1]) and NotRigid[w+1][1]=NotRigid[w][1] then
             NotRigid[w+1][2]:=NotRigid[w+1][2]-t;
             t:=t+1;
         else
             t:=1;
         fi;

    od;

    ##################################################################
    Boundary:=function(k,m)
        if m<0 then return NegateWord(BoundaryRec[k][AbsInt(m)]);fi;
        return BoundaryRec[k][m];
    end;

    Stabilizer:=function(k,m)
        return StabRec[k+1][m];
    end;

    Dimension:=function(k)
        if k>N then return 0;fi;
        return DimRec[k+1];
    end;

    Action:=function(k,i,j)
        return 1;
    end;
    ##################################################################
return Objectify(HapNonFreeResolution,
    rec(
    dimension:=Dimension,
    Partition:=Partition,
    boundary:=Boundary,
    homotopy:=fail,
    elts:=Elts,
    group:=C!.group,
    stabilizer:=Stabilizer,
    action:=Action,
    subdividing:=SubdividingCell,
    replacecell:=ReplaceCell,
    isrigid:=IsRigidCell,
    properties:=
    [["length",Maximum(1000,N)],
    ["characteristic",0],
    ["type","resolution"]]  ));
end);



################### end of ControlledSubdivision ############################

################## End of Hybrid Subdivision ##############################
DeclareGlobalFunction("VirtuallySimplicialSubdivision");
InstallGlobalFunction(VirtuallySimplicialSubdivision,
function(arg)
local W, StabRec, i, j, N, M, x, bdry, s1, s2, p, k, w, t,
      DimRec, BoundaryRec, id, dims, NotRigid, NewCell,
      Cell, Elts, Boundary, Dimension, CLeftCosetElt, LCoset,
      pos, Stab, Mult,DimTemp, BoundaryTemp, Partition, ConnectedCheck,
      Stabilizer, Action, IsRigidCell, ReplaceCell, SubdividingCell,
      Orbit, C, NaiveEulerChar, CellHomology, exportFile, IsContractible,
      CosetRec;

    C:=arg[1];
    i:=0;
    while C!.dimension(i)>0 do
        i:=i+1;
    od;
    N:=i-1; # Length of the chain complex
    M:=N;




if Length(arg)=2 then M:=Minimum(N,arg[2]);
fi;


    Elts:=C!.elts;
    StabRec:=[];
    DimRec:=[];
    CosetRec:=[];

    ##################################################################
    # If g in Elts return the position of g in the list,
    # otherwise, add g to Elts and return the position.
    pos:=function(g)
    local posit;

        posit:=Position(Elts,g);
        if posit=fail then
            Add(Elts,g);
            return Length(Elts);
        else
            return posit;
        fi;
    end;
    ##################################################################
    id:=pos(One(C!.group));
    ##################################################################
    # return the stabilizer of g*e,
    #
    Stab:=function(e,g)
    return ConjugateGroup(StabRec[e[1]+1][AbsInt(e[2])],Elts[g]^-1);
    end;
    ##################################################################
    # returns  a  "canonical"  representative  of  the  right  coset
    # Elts[g]*Stab[i+1][j]

    ######Add a record for left coset
    for i in [0..N] do
      CosetRec[i+1]:=[];
    od;

    CLeftCosetElt:=function(i,j,g)
      local p;
      if IsBound(CosetRec[AbsInt(i)+1][j]) then
        if IsBound(CosetRec[AbsInt(i)+1][j][g]) then
          return CosetRec[AbsInt(i)+1][j][g];
        else
          p:=pos(CanonicalRightCountableCosetElement
                                  (StabRec[AbsInt(i)+1][j],Elts[g]^-1)^-1);
          CosetRec[AbsInt(i)+1][j][g]:=p;
          return p;
        fi;
      else
        CosetRec[i+1][j]:=[];
        p:=pos(CanonicalRightCountableCosetElement
                                (StabRec[AbsInt(i)+1][j],Elts[g]^-1)^-1);
        CosetRec[AbsInt(i)+1][j][g]:=p;
        return p;
      fi;

    end;
    ##################################################################
    ##
    ##  Input:  A list L, degree k, position g of an element
    ##  Output: Product of g and L.
    ##
    Mult:=function(L,k,g)
    local x,w,t,h,y,vv,LL;
        vv:=[];
        LL:=ShallowCopy(L);
        for x in [1..Length(LL)] do
            w:=Elts[g]*Elts[LL[x][2]];
            t:=CLeftCosetElt(k,AbsInt(LL[x][1]),pos(w));
            Add(vv,[LL[x][1],t]);
        od;
        return vv;
    end;
    ###################################################################
    # Store essential data: stabilizers, boundaries, dimensions

    Orbit:=[];
    for i in [1..N+1] do
        Orbit[i]:=[];
    od;

    for i in [0..N] do
        StabRec[i+1]:=[];
        DimRec[i+1]:=C!.dimension(i);
        for j in [1..C!.dimension(i)] do
            StabRec[i+1][j]:=C!.stabilizer(i,j);
        od;
    od;

    BoundaryRec:=[];
    for i in [1..N] do
        BoundaryRec[i]:=[];
--> --------------------

--> maximum size reached

--> --------------------

[ zur Elbe Produktseite wechseln0.83Quellennavigators  Analyse erneut starten  ]