> # The aim of this little collection of prodecures is to provide the possibility of calculating the
> # homology groups of hypersurfaces of real toric varieties constructed with O. Viro's method of
> # combinatorial patchworking. These calculations can be done with the procedure "HS_Homology".
> # As a byproduct, I provide also a procedure "HS_NrOfComp" that counts the number of connected components only.
> # It has the advantage to be much faster than the homology calculations.
> # The Viro method requires (among others) a triangulation of a convex polytope as input data.
> # Obtaining such triangulations is tedious work if the dimension of the polytope and the number of simplices is not
> # small. I therefore include a tool that helps a little bit to manage this problem:
> # The procedures "prepare_simp" and "triangulate_simplex" compute a maximal triangulation of a simplex (the procedures were
> # written out of urgent need and are therefore neither elegant nor complete, as the user has to provide a first
> # (coarse) triangulation himself if he wants to use it with general polytopes).
> 
> # The procedures make use of the package "convex" provided by Matthias Franz, which can be found at
> # "http://www.math.uwo.ca/~mfranz/convex".
> 
> # Information on the mathematical background can be found in my thesis at "www.arlofin.de/Mathematik/Diss.html".
> # At the end of this file I include the examples which I worked on in the thesis and a couple of introductory examples
> # to illustrate the use of the procedures. Documentation of input parameters and output data
> # is also included in the comments of the source code. 
# Due to the need of runtime optimization I often make use of very compact Maple-typic expressions
# which are difficult to read especially for the inexperienced user. I am thinking on a "didactic"
> # version of this document, which would run slower but would make the code clearer...
> 
> # The procedures are free to any use, modification, distribution etc. If you found it useful I would be glad
> # to hear comments, bugs (hopefully not!), how you used it etc.
> 
> 
> # Florian Schwerteck, May 5 2010
> # Email: florian.schwerteck@arlofin.de
> 
> 
> 
> with(convex):
> 
> 
> ## ------------- Basic functions dealing with the homomorphisms "Lattice" --> {1, -1} --------------
> 
> ## The global variable LattHom2 is a list consisting of all d-dimensional +1/-1-vectors,
> ## representing the different homomorphisms in the canonical basis
## It must be defined prior to using this and the following procedures
> 
> Nr2Hom:=proc(n::integer)::list(integer);
> 
> ## Input: A number n
> ## Output: The n-th homomorphism stored in LattHom2
> global LattHom2;
> 
>    LattHom2[n]
> 
> end proc:
> 
> Hom2Nr:=proc(xi::list(integer))::integer;
> 
> ## Input: A homomorphism (a vector with +1/-1-entries)
> ## Output: A number wich indicates the position of the homomorphism in LattHom2
> ## The procedure is based upon the particular ordering of the homomomorphisms in LattHom2!!
> 
> option remember;
> 
> local number, i, add_hom, max_dim;
> 
> global LattHom2;
> 
> max_dim:=nops(LattHom2[1]);
> if nops(xi) <> max_dim then error("Dimension mismatch") end if;
> number:=1;
> number:=number+add(2^i*[0,1][xi[max_dim-i]], i=0..max_dim-1);  # Mult. Hom -> add. Hom.
> end proc:
> 
> LattHom2_mult:=proc(a::integer, b::integer)::integer;
> 
> ## Input: Numbers a and b
> ## Output: The number of the product (as homomorphisms) of the a-th with the b-th homomorphism in LattHom2
>  
>    option remember;
>    local a_hom, b_hom, i;
> 
>    global LattHom2;
> 
>    a_hom:= Nr2Hom(a);
>    b_hom:= Nr2Hom(b);
>    Hom2Nr([seq(a_hom[i]*b_hom[i], i=1..nops(a_hom))])
>    # In newer Maple versions equivalent to: Hom2Nr(a_hom *~ b_hom)
> end proc:
> 
> LattHom2_val:=proc(bas::list,expo::list)
> 
> 
> ## Input: A lattice element bas, which is assumed to be already reduced modulo 2
>           A homomorphism expo
> ## Output: The value of the homomorphism
>    option remember;
> 
>    local i, j, nr_of_signchanges;
> 
>    if nops(bas) <> nops(expo) then error("The lists must be of the same size") end if;
>    nr_of_signchanges:=nops(select(i -> not(bas[i] = 1 or irem(expo[i], 2) = 0), [seq(j, j=1..nops(bas))]));
>    if irem(nr_of_signchanges,2) <> 0 then -1
>       else 1
>    end if;
> end proc:
> 
> 
> Latt_basis:=proc(gen_Vecs_RE::{listlist(integer), list})::listlist(integer);
> 
> ##Input: A list of primitive lattice vectors in row echelon form, which span a Q-vector space V
> 
> ##Output: A list of lattice vectors (still in row echelon form) which form a basis of the lattice points in V
> 
>    uses LinearAlgebra:-Modular;
> 
>    local LBasis, V_dim, Amb_dim, i, lead_zeros, vec_nr, t_vec, g, primes, p, mult, factor_nr, Perm, Temp_Row, M, U, Q, rp, d, r, U2, U3, M_mod, t_vec_mod, M_mod_RE, k, t_vec_mod_RE, x, test, red_possible;
> 
>    if gen_Vecs_RE = [] then return [] end if;
>    if type(gen_Vecs_RE, listlist(integer)) = false then error "The argument should be a listlist(integer),if not empty" end if;
>    V_dim:=nops(gen_Vecs_RE);             ##Dimension of the subspace
>    Amb_dim:=nops(gen_Vecs_RE[1]);        ##Dimension of the ambient space
>    lead_zeros:=Array(1..V_dim);
>    for vec_nr from 1 to V_dim do
>       t_vec:=gen_Vecs_RE[vec_nr];
>       for i from 1 to Amb_dim while t_vec[i] = 0 do end do;
>       lead_zeros[vec_nr]:=i-1
>    end do;
>    LBasis:=Array(1..V_dim);
>    LBasis[V_dim]:=gen_Vecs_RE[V_dim];
>    for vec_nr from V_dim-1 to 1 by -1 do
>       t_vec:=gen_Vecs_RE[vec_nr];
>       M:=LinearAlgebra[Transpose](Matrix(V_dim-vec_nr, Amb_dim, convert(LBasis[vec_nr+1..V_dim], list)));
>       g:=igcd(op(t_vec[lead_zeros[vec_nr]+1..lead_zeros[vec_nr+1]]));
>       primes:=op(ifactors(g)[2..-1]);
>       for factor_nr from 1 to nops(primes) do
>          p, mult:=op(primes[factor_nr]);
>          M_mod:=Mod(p, M, integer[8]);
>          U:=Mod(p, M, integer[8]);               ##U:=M_mod would be horribly wrong!
>          t_vec_mod:=Mod(p, t_vec, integer[8]);
>          Q, rp, d:=RowEchelonTransform(p,U,true,true,true,false);   ##For the following see Maple Help on RET...
>          U:=LinearAlgebra[Transpose](U);
>          r := nops(rp);
>          U2 := Create(p, Amb_dim, Amb_dim, identity, integer[8]):
>          for i from 1 to r do U2[i]:=U[i] od:
>          U2:=LinearAlgebra[Transpose](U2);
>          Perm:=Create(p, Amb_dim, Amb_dim, identity, integer[8]):
>          for i from 1 to LinearAlgebra[Dimension](Q) do
>             Temp_Row:=Perm[i];
>             Perm[i]:=Perm[Q[i]];
>             Perm[Q[i]]:=Temp_Row
>          end do;
>          U3:=Multiply(p,U2, Perm);
>          M_mod_RE:=Multiply(p,U3, M_mod);                           ##...Here finally the RowEchelonForm
>          red_possible:=true;
>          for k from 1 to mult while red_possible do
>             t_vec_mod_RE:=Multiply(p,U3, t_vec_mod);
>             x:=Create(p,V_dim-vec_nr,1,integer[8]):
>             for i from 1 to V_dim-vec_nr do
>                x[rp[i]]:=t_vec_mod_RE[i]
>             end do;
>             test:=Multiply(p,M_mod_RE, x);
>             if evalb(convert(test[[1..-1],1], list) = convert(t_vec_mod_RE, list)) then         ##Sorry for these terrible but necessary type conversions! Complain to the Maple designers!
>                t_vec:=(t_vec-convert(M.x, list))/p;
>                t_vec_mod:=Mod(p, t_vec, integer[8]);
>             else
>                red_possible:=false
>             end if
>          end do
>       end do;
>       LBasis[vec_nr]:=t_vec
>    end do;
>    LBasis
> end proc:
> 
> Cosets:=proc(Subgroup::set(integer))::list;
> 
> ## Computes the cosets of a subgroup of LattHom2
> ## Input: The subgroup as a set of numbers indicating the positions in LattHom2
> ## Output: The cosets as list of sets of the above type
> 
> ## Attention: It is vital, but not verified by the procedure, that the input is indeed a subgroup!!
> 
> local N, n, unused, res,  Coset, Coset_nr, j;
> 
> global LattHom2;
> 
> N:=nops(LattHom2);
> n:=nops(Subgroup);
> if irem(N, n) <> 0 then error("Don't be stupid. This cannot be a subgroup") end if;
> res:=[Subgroup];
> unused:={seq(j,j=1..N)} minus Subgroup;
> for Coset_nr from 2 to N/n do
>    Coset:={seq(LattHom2_mult(Subgroup[j],op(1,unused)),j=1..n)};
>    unused:=unused minus Coset;
>    res:=[op(res),Coset]
> end do;   
> res
> end proc:
>  
> ## ----------- Basic functions dealing with matrices ----------------------
> 
> diagonalyze:=proc(M::{matrix, Matrix})
> 
> ## Input: A matrix
> ## Output: The number of ones on the diagonal, and the list of the remaining non-zero entries on the diagonal
>    local i, diag, diag_ones, eldiv, rdim, cdim, x;
> 
>    uses LinearAlgebra;
> 
>    rdim, cdim:=Dimension(M);
>    diag:=[seq(M[i,i], i=1..min(rdim, cdim))];
>    diag_ones, eldiv:=selectremove(x-> x=1, diag);
>    eldiv:=remove(x -> x=0, eldiv);
>    nops(diag_ones), eldiv
> end proc:
> 
> ## -----------Basic functions dealing with triangulations -----------------
> 
> calc_complx  := proc(l::{set, list})::set;
> 
> ## Input: A set or a list of (abstract) simplizes, which are in their turn
> ##        given as sets (non-degenerated) or lists (poss. degen.) of integers representing points
> ## Output: The set of all subsimplizes, given as sorted lists of points, with degen. removed 
> 
>    local simplices, simplex, i;
> 
>    simplices := {}:
>    for i from 1 to nops(l) do
>       if type(l[i], set) then
>          simplex:=l[i]
>       else
>          simplex:={op(l[i])}   #Removes degenerations
>       end if;
>       simplex:=sort(convert(simplex, list));
>       simplices:= simplices union {op(combinat[powerset](simplex))}     #Alg. preserves sorting! (at least with the currently implemented algorithms of powerset)
>    end do;
>    simplices minus {[]}
>    end proc:
> 
> 
> stratify := proc(Simplizes::{set, list}, dims::range)::Array;
> 
> ## Input: A set or a list of (abstract) simplizes, which are in their
> ##        turn given as sets (non-degenerated) or lists (poss. degen.) of points
> ## Output: An array, whose i-th entry is a list of the i-dimensional simplizes,
> ##        each of which is given as list of points.
> ##        If the second argument is given, then it determines the range of dimensions
> ##        of the output, otherwise, the maximal range is returned.
> 
> 
>    local Simplex_dimensions, min_dim, max_dim, Strata;
> 
>    
>    Simplex_dimensions:=seq(nops(Simplizes[i])-1, i=1..nops(Simplizes));
>    min_dim:=min(Simplex_dimensions);
>    max_dim:=max(Simplex_dimensions);
>    if _npassed = 2 then
>       min_dim:=max(op(1,_params[2]), min_dim);
>       max_dim:=min(op(2,_params[2]), max_dim);
>    end if;
>    Strata:=Array(min_dim..max_dim, [seq(convert(select(x -> nops(x)=i, Simplizes),list), i=min_dim+1..max_dim+1)])
>     ##The above expression may not be the fastest: It avoids a for-loop, but requires several passes over the data.
>     ##Tests necessary!
> end proc: 
> 
> checkTriangulation:=proc(A::listlist(integer), T::set)::boolean;
> 
> ## Checks, whether the simplizes in T define a triangulation of the
> ## convex hull of the points in A
> ## If not, the procedure produces an error.
> 
>    local P, dim_P, all_ok, T_max, Simplex, nr_Simplex, dim_Simplex, Simplex2, nr_Simplex2, dim_Simplex2, isSubSimplex, vol_T, vol_Simplex, IS;
> 
>    all_ok:=true;
>    P:=convhull(op(A));
>    dim_P:=dim(P);
>    if nops(A[1]) <> dim_P then
>       all_ok:=false;
>       error "The dimension of the polytope is too small";
>    end if;
> 
>    T_max:=[];
>    for Simplex in T do
>       dim_Simplex:= nops(Simplex)-1;
>       if dim_Simplex < dim_P then
>          isSubSimplex:=false;
>          for Simplex2 in T do
>             dim_Simplex2:=nops(Simplex2)-1;
>             if verify({op(Simplex)}, {op(Simplex2)}, `subset`) then
>                if dim_Simplex < dim_Simplex2 then
>                   isSubSimplex:=true
>                end if
>             end if
>          end do;
>          if isSubSimplex = false then
>             all_ok:=false;
>             error "%1 ist nicht in einem max. Simplex enthalten", `Simplex`;
>          end if
>       else
>          T_max:=[op(T_max), convhull(seq(A[i],i in Simplex))];
>          if dim(T_max[-1]) <> dim_Simplex then
>             all_ok :=false;
>             error "%1 ist kein richtiger Simplex", Simplex
>          end if
>       end if
>    end do;
> 
>    vol_T:=0;
>    for nr_Simplex from 1 to nops(T_max) do
>       Simplex:=T_max[nr_Simplex];
>       vol_Simplex:=volume(Simplex);
>       if vol_Simplex = 0 then
>          all_ok:=false;
>          error "%1 ist kein echter Simplex", `Simplex`;
>       end if;
>       vol_T:=vol_T + vol_Simplex;
>       for nr_Simplex2 from nr_Simplex +1 to nops(T_max) do
>          Simplex2:=T_max[nr_Simplex2];
>          IS:=intersection(Simplex, Simplex2);
>          if not (isface(IS, Simplex) and isface(IS, Simplex2)) then
>             all_ok:=false;
>             error "%1 und %2 schneiden sich nicht in einer gemeinsamen Seite", vertices(Simplex), vertices(Simplex2);
>          end if
>       end do
>    end do;
>    if vol_T <> volume(P) then
>       all_ok:=false;
>       error "Die Triangulierung überdeckt das Polytop nicht exakt"
>    end if;
>    all_ok
> end proc: 
> 
> 
> calc_ComplexIncidences:=proc(Complx::Array)::Array;
> 
> ## Input:  An array, whose d-th entry is a list of d-dimensional simplizes.
> ##         Simplizes are given as ordered list of numbers
> ## Output: An array, which describes the incidence relations between simplizes differing in their dimension by 1.
> ##         The d,n-th entry consists of two elements:
> ##           (1) The set of simplex numbers of simplizes in dimension
> ##               d+1, for which the n-th simplex in dimension d is a facet
> ##           (2) The list of simplex numbers of facets of the n-th simplex in dimension d,
> ##               such that the j-th entry in the list stores the information
> ##               on that facet obtained by deleting the j-th vertex
> 
>    option remember;
> 
>    local Complx_incidences, maxdim, mindim, sizeof_Complx, Simplex,
>          nr_Simplex, dim_Simplex, pSimplex,nr_pSimplex, j, parents;
> 
>    mindim:=max(op(1, ArrayDims(Complx)),0);
>    maxdim:=op(2, ArrayDims(Complx));
>    sizeof_Complx:=Array(mindim..maxdim, [seq(nops(Complx[j]), j=mindim..maxdim)]);
> 
>          ##We proceed by calculating the boundaries in dim+1, giving us also the "parents" in dim:
> 
>    Complx_incidences:=Array(mindim..maxdim);
>    Complx_incidences[mindim]:=Array(1..sizeof_Complx[mindim], [seq([], nr_Simplex=1..sizeof_Complx[mindim])]);     ##Minimal simps have no bd
>    
>    for dim_Simplex from mindim to maxdim-1 do
>       Complx_incidences[dim_Simplex+1]:=Array(1..sizeof_Complx[dim_Simplex+1], [seq(Array(1..dim_Simplex+2), i=1..sizeof_Complx[dim_Simplex+1])]);
>   
>       for nr_Simplex from 1 to sizeof_Complx[dim_Simplex] do
>          Simplex:=Complx[dim_Simplex][nr_Simplex];
>          parents:={};
>          for nr_pSimplex from 1 to sizeof_Complx[dim_Simplex+1] do
>             pSimplex:=Complx[dim_Simplex+1][nr_pSimplex];
>             for j from 1 to dim_Simplex+2 do
>                if Simplex = subsop(j = NULL, pSimplex) then
>                   parents:=parents union {nr_pSimplex};
>                   Complx_incidences[dim_Simplex+1][nr_pSimplex][j]:=nr_Simplex;
>                   break
>                end if
>             end do;
>          end do;
>          Complx_incidences[dim_Simplex][nr_Simplex]:= parents, convert(Complx_incidences[dim_Simplex][nr_Simplex], list)
>       end do
>    end do;
>    
>    for nr_Simplex from 1 to sizeof_Complx[maxdim] do
>       Complx_incidences[maxdim][nr_Simplex]:={},convert(Complx_incidences[maxdim][nr_Simplex], list)
>    end do;
>    Complx_incidences
> end proc:
> 
> ## -------------- Functions for the generation of a maximal triangulation (works with simplizes so far)
> 
> prepare_simp:=proc(A::listlist[integer])
> 
> local d, P, P_faces, P_faces_inc, A_strat, A_verts, A_faces, pt, Mf_dim, Mf_nr, face_dim, face_nr, face_verts, face_verts_pos, pos, vert_nr;
> 
> P:=convhull(op(A));
> d:=dim(P);
> if d <> nops(A[1]) then error "Sicherheitshalber bitte nur volldimensionale Simplizes" end if;
> if not issimplex(P) then error "Die Punkte müssen einen Simplex beschreiben!" end if;
> P_faces:=faces(P);
> P_faces_inc:=calc_FaceIncidences(P_faces);
> A_strat:=Array(0..d, [seq([],i=0..d)]);
> for pt in A do
>    Mf_dim,Mf_nr:=op(Minface([pt], P_faces, P_faces_inc, d,1));
>    A_strat[Mf_dim]:=[op(A_strat[Mf_dim]), [pt, Mf_nr]]
> end do;
> A_verts:=[seq(A_strat[0][i][1], i=1..nops(A_strat[0]))];
> A_faces:=Array(0..d, [seq([],i=0..d)]):
> for face_dim from 0 to d do
>    for face_nr from 1 to nops(P_faces[face_dim]) do
>       face_verts:=vertices(P_faces[face_dim][face_nr]);
>       face_verts_pos:=[];
>       for vert_nr from 1 to nops(face_verts) do
>          member(face_verts[vert_nr], A_verts,'pos');
>          face_verts_pos:=[op(face_verts_pos), pos]
>       end do;
>       A_faces[face_dim]:=[op(A_faces[face_dim]), [face_verts_pos, P_faces_inc[face_dim][face_nr]]]
>    end do:
> end do:
> 
> A_strat, A_faces
> 
> end proc:
> 
> triangulate_simplex:=proc(Simplex::Array, Faces::Array)::set;
> 
> local d,i,pt, pt_fnr, pt_fdim, pt_face_innerpts, nr, triang, face_verts, nr_faceverts, Subsimplex, ss_nr, k, del_pt_nr, fsimp_verts, A, P, L, U, R, Q, face_nrs, face_nrs_new, test_pt, v, x, x_neg, x_nonneg, zero_pos, plus_pos, j, testpt_fdim, face_found, found, test_face_nr, test_face;
> 
> d:=op(2, ArrayDims(Simplex));
> for i from d to 1 by -1 while nops(Simplex[i])=0 do end do;
> if i < 1 then
>    #if volume(convhull(seq(Simplex[0][i][1], i=1..nops(Simplex[0])))) > 1/d! then DEBUG() end if;
>    return Simplex
> else
>    pt_fdim:=i
> end if;
> pt, pt_fnr:=op(Simplex[pt_fdim][1]);
> pt_face_innerpts:={};
> for nr from 1 to nops(Simplex[pt_fdim]) do
>    if op(2, Simplex[pt_fdim][nr]) = pt_fnr then
>       pt_face_innerpts:=pt_face_innerpts union  {op(1, Simplex[pt_fdim][nr])}
>    end if
> end do;
> triang:={};
> face_verts:=Faces[pt_fdim][pt_fnr][1];
> nr_faceverts:=nops(face_verts);
> for ss_nr from 1 to nr_faceverts do
>    del_pt_nr:=Faces[pt_fdim][pt_fnr][1][ss_nr];
>    for k from 1 to nops(Faces[d-1]) while member(del_pt_nr, Faces[d-1][k][1]) do end do;
>    face_nrs:={k};
>    Subsimplex:=Array(0..d, [seq([], i=0..d)]);
>    for i from d-1 to 1 by -1 do
>       Subsimplex[i]:=select(x-> op(2,x) in face_nrs, Simplex[i]);
>       face_nrs_new:={};
>       for nr in face_nrs do
>          face_nrs_new:=face_nrs_new union Faces[i][nr][3]
>       end do;
>       face_nrs:=face_nrs_new;
>    end do;
>    Subsimplex[0]:=subsop(del_pt_nr=[pt, pt_fnr], Simplex[0]);
>    fsimp_verts:=subsop(ss_nr = NULL, face_verts);
>    A:=Matrix(nr_faceverts-1, d,[seq(Simplex[0][i][1]-pt, i in fsimp_verts)])^+;
>    P,L,U,R:=LinearAlgebra[LUDecomposition](A, method='RREF');
>    Q:=LinearAlgebra[MatrixInverse](P.L.U);
>    for test_pt in pt_face_innerpts minus {pt} do
>       v:=Vector[column](test_pt-pt);
>       x:=convert(Q.v, list)[1..pt_fdim];
>       x_neg, x_nonneg:=selectremove(t -> t < 0, x);
>       if nops(x_neg)=0 then
>          zero_pos, plus_pos:=selectremove(j -> x_nonneg[j] = 0, [seq(l, l=1..nops(x_nonneg))]);
>          for j from 1 to nops(zero_pos) do
>             if zero_pos[j] >=ss_nr then
>                zero_pos[j]:=zero_pos[j]+1
>             end if
>          end do;
>          plus_pos:=[op(select(x-> x < ss_nr, plus_pos)), ss_nr, op(select(x-> x >= ss_nr, plus_pos)+~1)];
>          test_face:=[seq(face_verts[i], i in plus_pos)];
>          testpt_fdim:=nops(test_face)-1;
>          #testpt_fdim:=pt_fdim-nops(zero_pos);
>          face_found:=false;
>          for test_face_nr from 1 to nops(Faces[testpt_fdim]) while face_found=false do
>             if test_face = Faces[testpt_fdim][test_face_nr][1] then # Kann sich Reihenfolge ändern?
>                face_found := true
>             end if  
>          end do;
>          if face_found = false then error "Kann nicht sein" end if;
>   #if not contains(convhull(seq(Subsimplex[0][i][1], i in test_face)), test_pt) then DEBUG() end if;      
>         Subsimplex[testpt_fdim]:=[op(Subsimplex[testpt_fdim]), [test_pt, test_face_nr -1]];
> #if not contains(convhull(seq(Simplex[0][i][1], i in face_verts)), test_pt) then DEBUG() end if;
>       end if
>    end do;
> 
>    triang:=triang union {Subsimplex}
> end do;
> op(map(triangulate_simplex, triang, Faces))
> 
> end proc:
> 
> ## -------------- Functions concerning polytopes
> 
> calc_FaceIncidences:=proc(Facearray::Array)::Array;
> 
> ## Input:  An array, whose d-th entry is a list of d-dimensional faces of a polyhedron.
> ##         
> ## Output: An array, which describes the incidence relations between faces differing in their dimension by 1.
> ##         The d,n-th entry consists of two sets:
> ##           (1) The set of face numbers of faces in dimension
> ##               d+1, for which the n-th face in dimension d is a facet
> ##           (2) The set of face numbers of facets of the n-th simplex in dimension d.
> 
>    local Face_incidences, j, maxdim, mindim, sizeof_Facearray, Face,
>          nr_Face, dim_Face, nr_Facet ,nr_testFace, Face_Facets, Facet_nos;
> 
>    mindim:=max(op(1, ArrayDims(Facearray)),0);
>    maxdim:=op(2, ArrayDims(Facearray));
>    sizeof_Facearray:=Array(mindim..maxdim, [seq(nops(Facearray[j]), j=mindim..maxdim)]);
>    
>    Face_incidences:=Array(mindim..maxdim);
>    Face_incidences[maxdim]:=Array(1..sizeof_Facearray[maxdim], [seq({}, nr_Face=1..sizeof_Facearray[maxdim])]);
> 
>    for dim_Face from maxdim to mindim+1 by -1 do
>       Face_incidences[dim_Face-1]:=Array(1..sizeof_Facearray[dim_Face-1], [seq({}, i=1..sizeof_Facearray[dim_Face-1])]);
>       
>       for nr_Face from 1 to sizeof_Facearray[dim_Face] do
>          Face:=Facearray[dim_Face][nr_Face];
>          Face_Facets:=pred(Face);
>          Facet_nos:={};
>          for nr_Facet from 1 to nops(Face_Facets) do  
>             for nr_testFace from 1 to sizeof_Facearray[dim_Face-1] do
>                if Face_Facets[nr_Facet] &= Facearray[dim_Face-1][nr_testFace] then
>                   Facet_nos:=Facet_nos union {nr_testFace};
>                   Face_incidences[dim_Face-1][nr_testFace]:=Face_incidences[dim_Face-1][nr_testFace] union {nr_Face}
>                end if
>             end do
>          end do;
>          Face_incidences[dim_Face][nr_Face]:= Face_incidences[dim_Face][nr_Face], Facet_nos
>       end do
>    end do;
>    
>    for nr_Face from 1 to sizeof_Facearray[mindim] do
>       Face_incidences[mindim][nr_Face]:=Face_incidences[mindim][nr_Face],{}
>    end do;
>    Face_incidences
> end proc:
> 
> Minface:=proc(Simplex::listlist(integer), Facearray::{array, Array}, Face_incidences::{array, Array}, start_dim::nonnegint, start_nr::posint)::list(integer);
> 
> ## Berechnet die Dimension und die Position der minimalen Seite von Facearray,
> ## die den Simplex enthält und echt in der Startseite enthalten ist. Falls es keine gibt,
> ## wird der Wert [start_dim, start_nr] zurückgegeben. Es wird also angenommen (wg. Performance
> ## aber nicht mehr überprüft!), dass
> ## der Simplex in der Startseite enthalten ist.
> ##   !!Proz. erfordert Paket convex!!
> 
>    local dim_range, Simp_MinFace, contained, max_dim, min_dim, testface_nos, face_dim, face_nr,  ambface_found, Face, verts_contained, i;
> 
> ## Beginn des Hauptteils:
>    
>    dim_range:=op(2, Facearray);
>    if nops([dim_range])>1 then
>       error("Das Seiten-Array hat zu viele Indexdimensionen")
>    end if;
>    min_dim, max_dim:=op(dim_range);
>    if start_dim > max_dim then error "Es gibt keine %1-dimensionalen Seiten", startdim end if;
>    Simp_MinFace:=[start_dim, start_nr];
>    contained:=true;
>    testface_nos:=Face_incidences[start_dim][start_nr][2];
>    for face_dim from start_dim-1 to min_dim by -1 while contained = true do
>       ambface_found:=false;
>       for face_nr in testface_nos while ambface_found = false do
>          Face:=convert(Facearray[face_dim][face_nr], POLYHEDRON);
>          verts_contained:=select(x-> contains(Face,x), Simplex);  ##avoids loop, but passes over complete data (seems ok, as most vertices should be contained in the face)
>          if nops(verts_contained) = nops(Simplex) then ambface_found:=true end if;
>          #verts_contained:=true: 
>          #for i from 1 to nops(Simplex) do
>             #verts_contained:=verts_contained and contains(Face, Simplex[i])
>          #end do;
>          #ambface_found:=verts_contained;
>          if ambface_found then
>             Simp_MinFace:=[face_dim, face_nr];
>             testface_nos:=Face_incidences[face_dim][face_nr][2]
>          end if;
>       end do;
>       contained:=ambface_found;
>    end do;
>    Simp_MinFace;         
> end proc:
> 
> calc_glueing:=proc(A::{listlist(integer)}, Triang::{array, Array}, Triang_inc::{array, Array})::Array;
> 
> ##     !!!Needs package "convex"!!!
> 
>  local triang_min_dim, triang_max_dim, dim_range, max_dim, face_dim, face_nr, Face, Face_verts, PFaces, P, PFaces_inc, Glue, ind_Glue, copy_ind, test_const, pface_nr, hom_nr, Face_base, Face_prebase, base_vec, i, step, Triang_MinFace, simp_dim, simp_nr, startdim, startnr, simp, simp_vertices, Triang_glue, Faceglueing, MinFace_dim, MinFace_nr;
> 
> global LattHom2;
> 
>  #-----------------Vorgeplänkel:
> 
>    dim_range:=ArrayDims(Triang);
>    triang_min_dim:=op(1, dim_range);
>    triang_max_dim:=op(2, dim_range);
>    if nops(A[1]) = 0 then error("Wasn das für ne Punktemenge???") end if;
>    P:=convhull(op(A));
>    max_dim:=dim(P);
>    if max_dim < nops(A[1]) then error("Die konvexe Hülle ist nicht volldimensional") end if;
>    PFaces:=faces(P);
>    PFaces_inc:=calc_FaceIncidences(PFaces);
> 
>    #---------------Minimale Seiten bestimmen:
> 
>         ##TMF[dim][i] ":=" [minface_dim, minface_nr] of the minimal face
>         ##in PFaces containing the i-th dim-dimensional simplex:
> 
>    Triang_MinFace:=Array(dim_range, [seq(Array(1..nops(Triang[d])) , d in dim_range)]);
>         ##Calculate the triang_max_dim case separately:
>    if triang_max_dim = max_dim then
>       Triang_MinFace[triang_max_dim]:=[seq([max_dim,1], i=1..nops(Triang[triang_max_dim]))]
>    else
>       for simp_nr from 1 to nops(Triang[triang_max_dim]) do
>          simp:=Triang[triang_max_dim][simp_nr];
>          if nops(simp) = 0 then
>             Triang_MinFace[triang_max_dim][simp_nr]:=[-1,1]   ## dealing with "[] <= []"
>          else
>             simp_vertices:=[seq(A[i], i in simp)];
>             Triang_MinFace[triang_max_dim][simp_nr]:=Minface(simp_vertices, PFaces, PFaces_inc, max_dim, 1);
>          end if;
>        end do;
>        Triang_MinFace[triang_max_dim]:=convert(Triang_MinFace[triang_max_dim], list)   ##For faster (reading) access
>    end if;
>    for simp_dim from triang_max_dim - 1 to triang_min_dim by -1 do
>       for simp_nr from 1 to nops(Triang[simp_dim]) do
>          simp:=Triang[simp_dim][simp_nr];
>          if nops(simp) = 0 then
>             Triang_MinFace[simp_dim][simp_nr]:=[-1,1]   ## dealing with "[] <= []"
>          else
>             simp_vertices:=[seq(A[i], i in simp)];
>             startdim, startnr:=max_dim, 1;
>             for i in Triang_inc[simp_dim][simp_nr][1] do
>                if op(1,Triang_MinFace[simp_dim+1][i]) < startdim then
>                   startdim, startnr:=op(Triang_MinFace[simp_dim+1][i])
>                end if
>             end do;
>             Triang_MinFace[simp_dim][simp_nr]:=Minface(simp_vertices, PFaces, PFaces_inc, startdim, startnr);
>          end if
>        end do;
>        Triang_MinFace[simp_dim]:=convert(Triang_MinFace[simp_dim], list)             ##For faster (reading) access
>    end do;
> 
>    #----------Verklebung der Seiten des Polytops bestimmen:
> 
>    Faceglueing:=Array(0..max_dim, [seq(Array(1..nops(PFaces[d])), d=0..max_dim)]);
>    for face_dim from max_dim to 0 by -1 do
>       for face_nr from 1 to nops(PFaces[face_dim]) do
>          Face:=PFaces[face_dim][face_nr];
>          ind_Glue:={1};
>          for pface_nr in PFaces_inc[face_dim][face_nr][1] do       ##Consider the glueing of the "parent" faces
>             for hom_nr in Faceglueing[face_dim+1][pface_nr] do
>                if not member(hom_nr, ind_Glue) then                ##Not yet glued?
>                   ind_Glue:=ind_Glue union map2(LattHom2_mult, hom_nr, ind_Glue) ##Glue it (i.e. the whole coset)!
>                end if
>             end do
>          end do;
>          Face_prebase:=lines(affinehull(convert(Face, POLYHEDRON)));
>          Face_base:=Latt_basis(Face_prebase);
>          Glue:=ind_Glue;
>          #Glue2:=ind_Glue;      
>          for copy_ind in {seq(i, i=1..nops(LattHom2))} minus ind_Glue do 
>             test_const := true;
>             for base_vec in Face_base while test_const = true do
>                if 1 <> LattHom2_val(LattHom2[copy_ind], base_vec) then test_const := false end if
>             end do;
>             if test_const = true then
>                Glue:=Glue union {copy_ind}
>             end if;
>          end do;
>          Faceglueing[face_dim][face_nr]:=Glue
>       end do;
>       Faceglueing[face_dim]:=convert(Faceglueing[face_dim], list)             ##For faster (reading) access
>    end do;
>    
>    #------------Verklebung auf die Simplizes anwenden:
> 
>    Triang_glue:=Array(dim_range, [seq(Array(1..nops(Triang[simp_dim]),0), simp_dim in dim_range)]);
>    # Triang_glue:=Array(dim_range, [seq([seq(Faceglueing[Triang_MinFace[simp_dim][simp_nr][1]][Triang_MinFace[simp_dim][simp_nr][2]], simp_dim=1..nops(Triang[simp_dim])], simp_dim in dim_range)]);
>    # Slower:
>    for simp_dim from triang_min_dim to triang_max_dim do
>       for simp_nr from 1 to nops(Triang[simp_dim]) do
>          MinFace_dim:=Triang_MinFace[simp_dim][simp_nr][1];
>          MinFace_nr:=Triang_MinFace[simp_dim][simp_nr][2];
>          Triang_glue[simp_dim][simp_nr]:=Faceglueing[MinFace_dim][MinFace_nr]
>       end do;
>       convert(Triang_glue[simp_dim], list)
>    end do;
>    Triang_glue;
> end proc:
> 
> ## ------------------- Functions concerning sign functions -------------------
> 
> isSimplexFull :=proc(Simplex::list(integer), vorz_fkt::list(integer))::boolean;
> 
> ##Berechnet, ob alle Ecken des Simplexes dasselbe Vorzeichen (-> false) haben oder nicht (->true).
> 
>    local i, equal, vorz;
> 
>    if nops(Simplex) < 2 then return false end if;
>    vorz:=vorz_fkt[Simplex[1]];
>    #equals:=select(x -> vorz = vorz_fkt[x], Simplex);   ## avoids loop, but passes over the complete data
>    #if nops(equals) = nops(Simplex) then false else true end if;
>    equal:=true;
>    for i from 2 to nops(Simplex) while equal=true do
>       equal:= evalb(vorz = vorz_fkt[Simplex[i]])
>    end do;
>    not equal;
> end proc:
> 
> SignFunc_expand:=proc(A::listlist(integer), Sign_Func::list(integer))::listlist(integer);
>    
>    local max_dim, nbr_pts, step, Sign_func_exp, copy_ind, i;
> 
>    global LattHom2;
> 
>    max_dim:=nops(A[1]);
>    nbr_pts:=nops(A);
>    if nops(Sign_Func) <> nbr_pts then error("Anzahl der Punkte und Vorzeichen verschieden") end if;
>    Sign_func_exp:=[seq([],i = 1..nops(LattHom2))];
>    for copy_ind from 1 to nops(LattHom2) do
>       for i from 1 to nbr_pts do                   # Sign_func_exp[copy_ind]:=Sign_Func*~[seq(LattHom2_val(LattHom2[copy_ind],A[i]), i=1..nbr_pts)] 
>          Sign_func_exp[copy_ind]:=[op(Sign_func_exp[copy_ind]),Sign_Func[i]*LattHom2_val(LattHom2[copy_ind],A[i])]
>       end do
>    end do;
>    Sign_func_exp;
> end proc:
> 
> ## ------------------- Other functions --------------------------------
> 
> calc_OrientStruct:=proc(Sign_func::list(integer))::list;
> 
>    local v0,  gamma, plus_nr, minus_nr, sign_sep, alpha, i;
> 
>    plus_nr:=0;
>    minus_nr:=0;
>    sign_sep:=Array(1..nops(Sign_func));
>    
>    v0:=Sign_func[1];
>    for i from 1 to nops(Sign_func) do
>       if Sign_func[i] = v0 then
>          plus_nr:=plus_nr+1;
>          sign_sep[i]:=plus_nr;
>       else
>          minus_nr:=minus_nr+1;
>          sign_sep[i]:=-minus_nr;
>       end if
>    end do;
>    [plus_nr-1, minus_nr-1, convert(sign_sep, list)]
> end proc:
> 
> ## ------------ Main functions --------------------------------
> 
> HS_Homology:=proc(A::listlist(integer), T::set, Sign_Func, {dim_range::range := 0..infinity, check_T::boolean:=true, Euler::boolean:=false, coeff::name:=integer})::Array;
> 
>    local modulus, Triang, Triang_incidences, dim_P, min_dim, max_dim, Glueing_structure, sign_nr, Signs, Sign_Func_proc, Signs_global, HS_Cells, HS_Cells_StartIndex, HS_Cells_OrientStructure, HS_Cells_BoundaryMaps, simp_dim, simp_nr, cell_dim, cell_nr, Copies, Coset, i,j, k, bm_min_dim, bm_max_dim, bd_pos, Cell, bsimp_nr, bcell_nr, bCell, boundary_found, sign_bd_pos, alpha, beta, Vert_signs, corr, Homolgy, h_dim, rdim, cdim, Smith_Matrix, nb_of_ones, eldiv, rk, results;
> 
>    global LattHom2;
> 
>    if coeff=integer then
>       modulus := 0
>    elif StringTools[IsPrefix]("mod", convert(coeff, string)) then
>       modulus:=convert(coeff, string)[4..-1];
>       if StringTools[IsDigit](modulus) then
>          modulus:=parse(modulus);
>          if type(modulus, prime) = false or modulus > 127 then
>             error "The modulus of the coeffficient must be a prime <= 127"
>          end if
>       else
>          error "The modulus of the coeffficients must be a number (even a prime) and not %1%", modulus
>       end if
>    else
>       error "Unknown coefficients (must be integer or modp)"
>    end if;
>    
>    dim_P:=nops(A[1]);
>    if dim_P < 2 then error("Das Problem hat zu wenig Dimension") end if;
>    min_dim:=max(0, op(1, dim_range)-1);
>    max_dim:=min(dim_P -1, op(2, dim_range)+1);
>    if Euler = true then
>       min_dim:=0;
>       max_dim:=dim_P -1
>    end if;
>    
>    if check_T = true then
>       print("Checking the triangulation");
>       checkTriangulation(A, T)
>    end if;
>    
>    Triang := stratify(calc_complx(T)) [max(1, min_dim)..min(dim_P, max_dim+1)];  # dim_range in stratify benutzen!
> 
>    print("Calculating the incidences of the triangulation");
>    Triang_incidences:=calc_ComplexIncidences(Triang);
> 
>    print("Calculating the glueing structure");
>    LattHom2:=[[]];
>    for i from 1 to dim_P do
>       LattHom2:=[seq([1,op(k)], k in LattHom2), seq([-1,op(k)], k in LattHom2)]
>    end do;
>    Glueing_structure:=calc_glueing(A, Triang, Triang_incidences);
> 
>    if type(Sign_Func, list(integer)) then
>       Sign_Func_proc:=proc(i)
>          if i=1 then Sign_Func
>          else [] end if
>       end proc
>    elif type(Sign_Func, procedure) then
>       Sign_Func_proc:=Sign_Func;
>       print("Calculating the rest")
>    else
>       error "Sign_Func must either be a list of +/- 1 or a procedure generating such lists"
>    end if;
>    sign_nr:=1;
>    Signs:=Sign_Func_proc(1);
>    results:={};
>    while Signs <> [] do
> 
> # ---------Cells and orientations -----------
> 
>       if type(Sign_Func, list(integer)) then
>          print("Determine the cells of the hypersurface and their orientation")
>       end if;
> 
>       Signs_global:=SignFunc_expand(A,Signs);
>       HS_Cells:=Array(min_dim..max_dim);
>       HS_Cells_StartIndex:=Array(min_dim..max_dim);
>       HS_Cells_OrientStructure:=Array(min_dim..max_dim);
>       for cell_dim from min_dim to max_dim do
>          HS_Cells[cell_dim]:=[];
>          HS_Cells_StartIndex[cell_dim]:=[];
>          HS_Cells_OrientStructure[cell_dim]:=[];
>          cell_nr:=1;
>          for simp_nr from 1 to nops(Triang[cell_dim+1]) do
>             HS_Cells_StartIndex[cell_dim]:=[op(HS_Cells_StartIndex[cell_dim]), cell_nr];
>             Copies:=Cosets(Glueing_structure[cell_dim+1][simp_nr]);
>             for Coset in Copies do
>                if isSimplexFull(Triang[cell_dim+1][simp_nr], Signs_global[op(1,Coset)]) then
>                   HS_Cells[cell_dim]:=[op(HS_Cells[cell_dim]), [simp_nr, Coset]];
>                   HS_Cells_OrientStructure[cell_dim]:=[op(HS_Cells_OrientStructure[cell_dim]), calc_OrientStruct([seq(Signs_global[op(1, Coset)][i], i in Triang[cell_dim+1][simp_nr])])];
>                   cell_nr:=cell_nr+1
>                end if
>             end do
>          end do
>       end do;
> 
>       # -----------Boundary maps -----------------------------
> 
>       if type(Sign_Func, list(integer)) then
>          print("Calculating the boundary maps")
>       end if;
> 
>       bm_min_dim:=max(1,min_dim);
>       bm_max_dim:=max_dim;
>       HS_Cells_BoundaryMaps:=Array(bm_min_dim..bm_max_dim);
>       for cell_dim from bm_min_dim to bm_max_dim do
>          rdim:=nops(HS_Cells[cell_dim]);
>          cdim:=nops(HS_Cells[cell_dim-1]);
>          if modulus = 0 then
>             HS_Cells_BoundaryMaps[cell_dim]:=Matrix(rdim, cdim, storage=sparse)
>          else
>             HS_Cells_BoundaryMaps[cell_dim]:=LinearAlgebra:-Modular[Create](modulus, rdim, cdim, 0, integer[1])
>          end if;
>          for cell_nr from 1 to rdim do
>             Cell:=HS_Cells[cell_dim][cell_nr];
>             simp_nr:=Cell[1];
>             alpha, beta, Vert_signs:=op(HS_Cells_OrientStructure[cell_dim][cell_nr]);
>             for bd_pos from 1 to nops(Triang_incidences[cell_dim+1][simp_nr][2]) do
>                boundary_found:=false;
>                bsimp_nr:=Triang_incidences[cell_dim+1][simp_nr][2][bd_pos];
>                bcell_nr:=HS_Cells_StartIndex[cell_dim-1][bsimp_nr];
>                if bcell_nr > nops(HS_Cells[cell_dim-1]) then next end if;
>                bCell:=HS_Cells[cell_dim-1][bcell_nr];
>                while boundary_found = false and bCell[1] = bsimp_nr do
>                   boundary_found:=verify(Cell[2], bCell[2], `subset`);
>                   if not boundary_found then
>                      bcell_nr:=bcell_nr+1;
>                      if bcell_nr > nops(HS_Cells[cell_dim-1]) then break end if;
>                      bCell:=HS_Cells[cell_dim-1][bcell_nr]
>                   end if
>                end do;
>                if boundary_found then
>                   sign_bd_pos:=Vert_signs[bd_pos]-1;
>                   if bd_pos=1 and Vert_signs[2] < 0 then
>                      corr:=(alpha-1)*beta
>                   else
>                      corr:=0
>                   end if;
>                   if sign_bd_pos < 0 then
>                      if irem(sign_bd_pos,2) = 0 then
>                         HS_Cells_BoundaryMaps[cell_dim][cell_nr,bcell_nr]:=1
>                      else
>                         HS_Cells_BoundaryMaps[cell_dim][cell_nr,bcell_nr]:=-1
>                      end if
>                   else
>                      if irem(beta + corr + sign_bd_pos, 2) = 0 then
>                         HS_Cells_BoundaryMaps[cell_dim][cell_nr,bcell_nr]:=1
>                      else
>                         HS_Cells_BoundaryMaps[cell_dim][cell_nr,bcell_nr]:=-1
>                      end if
>                   end if
>                end if
>             end do
>          end do
>       end do;
>   
>       # ------------The Smith normal forms------------------
> 
>       if type(Sign_Func, list(integer)) then
>          if modulus=0 then
>             print("Calculating the Smith normal forms")
>          else
>             print("Calculating the ranks")
>          end if
>       end if;
> 
>       min_dim:=max(0, op(1, dim_range));
>       max_dim:=min(dim_P -1, op(2, dim_range));
>       Homolgy:=Array(min_dim..max_dim);
>       Homolgy[min_dim]:=nops(HS_Cells[min_dim]);
>       if modulus = 0 then
>          if min_dim > 0 then
>             Homolgy[min_dim]:=Homolgy[min_dim]-LinearAlgebra[Rank](HS_Cells_BoundaryMaps[min_dim])
>          end if;  # min_dim=0 -> Rank(..) = 0
>          for h_dim from min_dim+1 to max_dim do
>             Smith_Matrix:=LinearAlgebra[SmithForm](HS_Cells_BoundaryMaps[h_dim], method=integer);
>             nb_of_ones, eldiv:=diagonalyze(Smith_Matrix);
>             rk:=nb_of_ones+nops(eldiv);
>             Homolgy[h_dim]:=nops(HS_Cells[h_dim])-rk;
>             Homolgy[h_dim -1]:=modz(Homolgy[h_dim -1]-rk,eldiv)
>          end do;
>          if max_dim = dim_P-1 then
>             Homolgy[max_dim]:=modz(Homolgy[max_dim],[])
>          elif max_dim = 0 then
>             rk:=LinearAlgebra[Rank](HS_Cells_BoundaryMaps[max_dim+1]);
>             Homolgy[max_dim]:=modz(Homolgy[max_dim]-rk,[])
>          else
>             Smith_Matrix:=LinearAlgebra[SmithForm](HS_Cells_BoundaryMaps[max_dim+1], method=integer);
>             nb_of_ones, eldiv:=diagonalyze(Smith_Matrix);
>             rk:=nb_of_ones+nops(eldiv);
>             Homolgy[max_dim]:=modz(Homolgy[max_dim]-rk,eldiv)
>          end
>       else
>          if min_dim > 0 then
>             Homolgy[min_dim]:=Homolgy[min_dim]-LinearAlgebra:-Modular[Rank](modulus, HS_Cells_BoundaryMaps[min_dim], REF)
>          end if;  # min_dim=0 -> Rank(..) = 0
>          for h_dim from min_dim+1 to max_dim do
>             rk:=LinearAlgebra:-Modular[Rank](modulus, HS_Cells_BoundaryMaps[h_dim], REF);
>             Homolgy[h_dim]:=nops(HS_Cells[h_dim])-rk;
>             Homolgy[h_dim -1]:=modz(0, [seq(modulus, i=1..Homolgy[h_dim -1]-rk)])
>          end do;
>          if max_dim = dim_P-1 then
>             Homolgy[max_dim]:=modz(0, [seq(modulus, i=1..Homolgy[max_dim])])
>          else
>             rk:=LinearAlgebra:-Modular[Rank](modulus, HS_Cells_BoundaryMaps[max_dim+1], REF);
>             Homolgy[max_dim]:=modz(0, [seq(modulus, i=1..Homolgy[max_dim]-rk)])
>          end if
>       end if;
>       if Euler = true then
>          results:=results union {[chi=add((-1)^i*nops(HS_Cells[i]), i=0..dim_P-1), seq(h||h_dim = Homolgy[h_dim], h_dim = min_dim .. max_dim)]}
>       else
>          results:=results union {[seq(h||h_dim = Homolgy[h_dim], h_dim = min_dim .. max_dim)]}
>       end if;
>       sign_nr:=sign_nr+1;
>       Signs:=Sign_Func_proc(sign_nr)
>    end do;      ##->repeat with next sign function
> 
> 
>    if type(Sign_Func, list) then
>       op(results)
>    else
>       results
>    end if
> end proc:
> 
> 
> HS_NrOfComp:=proc(A::listlist(integer), T::set, Sign_Func, {smooth::boolean:=false, check_T::boolean:=true, Euler::boolean:=false})::nonnegint;
> 
>    local Triang, Triang_incidences, dim_P, min_dim, max_dim, Glueing_structure,  Sign_Func_proc, sign_nr, Signs, Signs_global, results, HS_Cells, HS_Cells_StartIndex, HS_Cells_BoundaryMaps, simp_dim, simp_nr, cell_dim, cell_nr, Copies, Coset, i,j, bd_pos, Cell, bsimp_nr, bcell_nr, bCell, boundaries, boundary_found, Components, Bd_Comp_Intersections, Intersection_indices, NoIntersection_indices, nrof_Intersections, comp_nr;
> 
>    global LattHom2;
> 
>    dim_P:=nops(A[1]);
>    if dim_P < 2 then error("Das Problem hat zu wenig Dimension") end if;
>    if Euler = true then
>       min_dim:=0;
>       max_dim:=dim_P -1
>    elif smooth = true then
>       min_dim:=dim_P-2;
>       max_dim:=dim_P-1
>    else
>       min_dim:=0;
>       max_dim:=1
>    end if;
>    
>    if check_T = true then
>       print("Checking the triangulation");
>       checkTriangulation(A, T)
>    end if;
> 
>    Triang := stratify(calc_complx(T))[min_dim+1..max_dim+1];
> 
>    print("Calculating the incidences of the triangulation");
>    Triang_incidences:=calc_ComplexIncidences(Triang);
> 
>    print("Calculating the glueing structure");
>    LattHom2:=[[]];
>    for i from 1 to dim_P do
>       LattHom2:=[seq([1,op(k)], k in LattHom2), seq([-1,op(k)], k in LattHom2)]
>    end do;
>    Glueing_structure:=calc_glueing(A, Triang, Triang_incidences);
> 
>    if type(Sign_Func, list(integer)) then
>       Sign_Func_proc:=proc(i)
>          if i=1 then Sign_Func
>          else [] end if
>       end proc
>    elif type(Sign_Func, procedure) then
>       Sign_Func_proc:=Sign_Func;
>       print("Calculating the rest")
>    else
>       error "Sign_Func must either be a list of +/- 1 or a procedure generating such lists"
>    end if;
>    sign_nr:=1;
>    Signs:=Sign_Func_proc(1);
>    results:={};
>    while Signs <> [] do
>    
>       if type(Sign_Func, list(integer)) then
>          print("Determine the cells of the hypersurface")
>       end if;
> 
>       Signs_global:=SignFunc_expand(A,Signs);
>       HS_Cells:=Array(min_dim..max_dim);
>       HS_Cells_StartIndex:=Array(min_dim..max_dim);
>       for cell_dim from min_dim to max_dim do
>          HS_Cells[cell_dim]:=[];
>          HS_Cells_StartIndex[cell_dim]:=[];
>          cell_nr:=1;
>          for simp_nr from 1 to nops(Triang[cell_dim+1]) do
>             HS_Cells_StartIndex[cell_dim]:=[op(HS_Cells_StartIndex[cell_dim]), cell_nr];
>             Copies:=Cosets(Glueing_structure[cell_dim+1][simp_nr]);
>             for Coset in Copies do
>                if isSimplexFull(Triang[cell_dim+1][simp_nr], Signs_global[op(1,Coset)]) then
>                   HS_Cells[cell_dim]:=[op(HS_Cells[cell_dim]), [simp_nr, Coset]];
>                    cell_nr:=cell_nr+1
>                end if
>             end do
>          end do
>       end do;
> 
>       if type(Sign_Func, list(integer)) then
>          print("Calculating the number of components")
>       end if;
>       
>       
>       if smooth then
>          cell_dim :=dim_P-1
>       else
>          cell_dim:=1
>       end if;   
>       Components:=[];     ##Each Component will be given as a set of its 0-cells (codim 1-cells in the smooth case)
>       for cell_nr from 1 to nops(HS_Cells[cell_dim]) do
>          boundaries:={};
>          Cell:=HS_Cells[cell_dim][cell_nr];
>          simp_nr:=Cell[1];
>          for bd_pos from 1 to nops(Triang_incidences[cell_dim+1][simp_nr][2]) do
>             boundary_found:=false;
>             bsimp_nr:=Triang_incidences[cell_dim+1][simp_nr][2][bd_pos];
>             bcell_nr:=HS_Cells_StartIndex[cell_dim-1][bsimp_nr];
>             if bcell_nr > nops(HS_Cells[cell_dim-1]) then next end if;
>             bCell:=HS_Cells[cell_dim-1][bcell_nr];
>             while boundary_found = false and bCell[1] = bsimp_nr do
>                boundary_found:=verify(Cell[2], bCell[2], `subset`);
>                if not boundary_found then
>                   bcell_nr:=bcell_nr+1;
>                   if bcell_nr > nops(HS_Cells[cell_dim-1]) then break end if;
>                   bCell:=HS_Cells[cell_dim-1][bcell_nr]
>                end if
>             end do;
>             if boundary_found then
>                boundaries:=boundaries union {bcell_nr}
>             end if
>          end do;
> 
>          Bd_Comp_Intersections:=map(`intersect`, Components, boundaries);
>          Intersection_indices, NoIntersection_indices:=selectremove(j-> Bd_Comp_Intersections[j]<>{}, [seq(i, i=1..nops(Components))]);
>          nrof_Intersections:=nops(Intersection_indices);
> 
>          if nrof_Intersections = 0 then
>             Components:=[op(Components), boundaries]
>          else
>             Components:=[`union`(boundaries, seq(Components[Intersection_indices[i]], i=1..nrof_Intersections)), seq(Components[NoIntersection_indices[i]], i=1..nops(NoIntersection_indices))]
>          end if
>       end do;
>         
>       if Euler = true then
>          results:=results union {[chi=add((-1)^i*nops(HS_Cells[i]), i=0..dim_P-1), Nr=nops(Components)]}
>       else
>          results:=results union {nops(Components)}
>       end if;
>       sign_nr:=sign_nr+1;
>       Signs:=Sign_Func_proc(sign_nr)
>    end do;
>    
>    if type(Sign_Func, list) then
>       op(results)
>    else
>       results
>    end if
> end proc:
> 
> 
> 
> with(convex):
Convex version 1.1.3, Copyright (C) 1999-2009 Matthias Franz
This package is distributed under the GNU General Public License,
see http://www-fourier.ujf-grenoble.fr/~franz/convex/ for more information.
> A := [[0, 0], [0, 1], [1, 0], [1, 1]];
                      [[0, 0], [0, 1], [1, 0], [1, 1]]
> T := {{1, 2, 4} , {4, 3, 1}};
                           {{1, 2, 4}, {1, 3, 4}}
> epsilon:=[1,-1,1,1];
                                [1, -1, 1, 1]
> h:=HS_Homology(A, T, epsilon, check_T = false, coeff=integer);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
              [h0 = _OBJ(MODZ, 1, []), h1 = _OBJ(MODZ, 1, [])]
> stopat(HS_Homology, 57);
                                [HS_Homology]
> unstopat();
                                     []
> nC:=HS_NrOfComp(A, T, epsilon, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      1
> epi:=proc(i)
> local vorz, j,k, a2m, res, anz_pkt;
> 
> anz_pkt:=4;
> a2m:=[1,-1];
> res:=[];
> if i <= 2^anz_pkt then
>    j:=i;
>    for k from 1 to anz_pkt do
>    vorz:=a2m[irem(j, 2, 'j')+1];
>    res:=[op(res), vorz];
>    end do
> end if;
> res
> end proc:
> gen_epi:=proc(anz_pkt::posint)
> 
> global epi;
> 
> epi:=proc(i)
> local vorz, j,k, a2m, res;
> 
> a2m:=[1,-1];
> res:=[];
> if i <= 2^anz_pkt then
>    j:=i;
>    for k from 1 to anz_pkt do
>    vorz:=a2m[irem(j, 2, 'j')+1];
>    res:=[op(res), vorz];
>    end do
> end if;
> res
> end proc;
> 
> end proc:
> gen_epi(4);
> epi(5);
proc(i)  ...  end;
                               [-1, 1, -1, 1]
> hh:=HS_Homology(A, T, epi, check_T = false, Euler = true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                           "Calculating the rest"
         {[chi = 0, h0 = _OBJ(MODZ, 1, []), h1 = _OBJ(MODZ, 1, [])]}
> nnrC:=HS_NrOfComp(A, T, epi, check_T = false, Euler = true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                           "Calculating the rest"
                             {[chi = 0, Nr = 1]}
> B := [[0, 0], [0, 2], [2, 0], [2, 2]];
                      [[0, 0], [0, 2], [2, 0], [2, 2]]
> 
> T2:={{1,2,3},{2,3,4}};
                           {{1, 2, 3}, {2, 3, 4}}
> eps2:=[-1,1,1,-1];
                               [-1, 1, 1, -1]
> h:=HS_Homology(B, T2, eps2);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
              [h0 = _OBJ(MODZ, 2, []), h1 = _OBJ(MODZ, 2, [])]
> nC:=HS_NrOfComp(B, T2, eps2);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      2
> 
> B_surf:=[seq([op(B[i]),0],i=1..nops(B)), [0,0,2]];
           [[0, 0, 0], [0, 2, 0], [2, 0, 0], [2, 2, 0], [0, 0, 2]]
> eB_surfp:=[op(eps2), 1];
                              [-1, 1, 1, -1, 1]
> TB_surf:={seq(T2[i] union {5}, i=1..nops(T2))};
                        {{1, 2, 3, 5}, {2, 3, 4, 5}}
> hB_surf:=HS_Homology(B_surf, TB_surf, eB_surfp, check_T=false, Euler=true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
         [chi = 4, h0 = _OBJ(MODZ, 2, []), h1 = _OBJ(MODZ, 0, []), 

           h2 = _OBJ(MODZ, 2, [])]
> 
> 
> Q:=[];
> for i from 0 to 4 do
>    for j from 0 to 4 do
>       Q:=[op(Q), [i,j]]
>    end do
> end do:
                                     []
> Q;
 [[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [1, 0], [1, 1], [1, 2], [1, 3], 

   [1, 4], [2, 0], [2, 1], [2, 2], [2, 3], [2, 4], [3, 0], [3, 1], [3, 2], 

   [3, 3], [3, 4], [4, 0], [4, 1], [4, 2], [4, 3], [4, 4]]
> TQ1:={[1,2,10],[2,3,10],[3,4,10],[4,5,10],[1,9,10],[1,8,9],[1,6,7],[1,8,14],[1,14,20],[1,7,20],[7,13,20],[6,7,13],[6,13,19],[6,19,25],[6,12,25],[6,11,12],[11,12,18],[11,17,18],[11,16,17],[12,18,25],[13,19,20],[14,15,20],[8,14,15],[8,9,15],[9,10,15],[16,21,22],[16,22,23],[16,23,24],[16,24,25],[16,17,25],[17,18,25],[19,20,25]};
{[1, 2, 10], [1, 6, 7], [1, 7, 20], [1, 8, 9], [1, 8, 14], [1, 9, 10], 

  [1, 14, 20], [2, 3, 10], [3, 4, 10], [4, 5, 10], [6, 7, 13], [6, 11, 12], 

  [6, 12, 25], [6, 13, 19], [6, 19, 25], [7, 13, 20], [8, 9, 15], [8, 14, 15], 

  [9, 10, 15], [11, 12, 18], [11, 16, 17], [11, 17, 18], [12, 18, 25], 

  [13, 19, 20], [14, 15, 20], [16, 17, 25], [16, 21, 22], [16, 22, 23], 

  [16, 23, 24], [16, 24, 25], [17, 18, 25], [19, 20, 25]}
> eQ1:=[1,1,1,1,1,1,-1,1,-1,1,1,1,-1,1,1,1,-1,1,-1,1,1,1,1,1,1];
[1, 1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 1

  ]
> hQ1:=HS_Homology(Q, TQ1, eQ1);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
             [h0 = _OBJ(MODZ, 10, []), h1 = _OBJ(MODZ, 10, [])]
> nrCQ1:=HS_NrOfComp(Q, TQ1, eQ1, check_T=false); 
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                     10
> nops(Q);
                                     25
> gen_epi(25):
> nC_multQ1:=HS_NrOfComp(Q, TQ1, epi, check_T=false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                           "Calculating the rest"
Warning,  computation interrupted
> P:=[];
> for i from 0 to 6 do
>    for j from 0 to 6-i do
>       P:=[op(P), [i,j]]
>    end do
> end do:
                                     []
> P;
 [[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [0, 5], [0, 6], [1, 0], [1, 1], 

   [1, 2], [1, 3], [1, 4], [1, 5], [2, 0], [2, 1], [2, 2], [2, 3], [2, 4], 

   [3, 0], [3, 1], [3, 2], [3, 3], [4, 0], [4, 1], [4, 2], [5, 0], [5, 1], 

   [6, 0]]
> TP1:={[1,2,8],[2,8,9],[2,3,9],[3,9,10],[3,4,10],[4,5,10],[5,10,15],[5,11,15],[5,6,11],[6,11,19],[6,16,19],[6,12,16],[6,7,12],[7,12,13],[8,9,14],[10,14,15],[9,10,14],[11,15,19],[12,16,17],[12,13,17],[13,17,18],[14,15,19],[16,19,20],[16,17,20],[17,20,21],[17,18,21],[18,21,22],[19,20,23],[20,21,23],[21,23,24],[21,24,25],[21,22,25],[23,24,26],[24,25,26],[25,26,27],[26,27,28]};
{[1, 2, 8], [2, 3, 9], [2, 8, 9], [3, 4, 10], [3, 9, 10], [4, 5, 10], 

  [5, 6, 11], [5, 10, 15], [5, 11, 15], [6, 7, 12], [6, 11, 19], [6, 12, 16], 

  [6, 16, 19], [7, 12, 13], [8, 9, 14], [9, 10, 14], [10, 14, 15], 

  [11, 15, 19], [12, 13, 17], [12, 16, 17], [13, 17, 18], [14, 15, 19], 

  [16, 17, 20], [16, 19, 20], [17, 18, 21], [17, 20, 21], [18, 21, 22], 

  [19, 20, 23], [20, 21, 23], [21, 22, 25], [21, 23, 24], [21, 24, 25], 

  [23, 24, 26], [24, 25, 26], [25, 26, 27], [26, 27, 28]}
> eP1:=[1,1,1,1,1,1,-1,1,-1,1,-1,1,1,1,1,-1,1,-1,1,1,1,1,-1,-1,-1,1,1,-1];
[1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 

  -1, -1, 1, 1, -1]
> hP1:=HS_Homology(P, TP1,eP1);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
             [h0 = _OBJ(MODZ, 11, []), h1 = _OBJ(MODZ, 11, [])]
> hhP1:=HS_NrOfComp(P, TP1, eP1, check_T=false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                     11

> V:=[[-2,3],[-1,4],[0,4],[1,4],[-1,5],[0,5],[1,5],[0,6],[1,6],[1,7],[2,7],[2,8],[3,9]];
[[-2, 3], [-1, 4], [0, 4], [1, 4], [-1, 5], [0, 5], [1, 5], [0, 6], [1, 6], 

  [1, 7], [2, 7], [2, 8], [3, 9]]
> TV1:={[1,2, 5],[1,2,3],[1,3,4],[2,5,8],[2,6,8],[2,3,6],[3,4,6],[4,6,7],[4,7,11],[4,11,13],[6,8,10],[6,9,10],[6,7,9],[7,9,11],[9,10,11],[10,11,12], [11,12,13]};
{[1, 2, 3], [1, 2, 5], [1, 3, 4], [2, 3, 6], [2, 5, 8], [2, 6, 8], [3, 4, 6], 

  [4, 6, 7], [4, 7, 11], [4, 11, 13], [6, 7, 9], [6, 8, 10], [6, 9, 10], 

  [7, 9, 11], [9, 10, 11], [10, 11, 12], [11, 12, 13]}
> eV1:=[1,-1,1,-1,1,1,-1,-1,1,1,-1,1,1];
                [1, -1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1]
> hV:=HS_Homology(V, TV1, eV1, check_T = false, Euler = true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
          [chi = 0, h0 = _OBJ(MODZ, 4, []), h1 = _OBJ(MODZ, 4, [])]
> hhV:=HS_NrOfComp(V, TV1, eV1, check_T=false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      4
> 
> Q_surf:=[]:
> for pkt_nr from 1 to nops(Q) do
>    Q_surf:=[op(Q_surf), [op(Q[pkt_nr]),0]]
> end do:
> Q_surf:=[op(Q_surf), [0,0,2]];
[[0, 0, 0], [0, 1, 0], [0, 2, 0], [0, 3, 0], [0, 4, 0], [1, 0, 0], [1, 1, 0], 

  [1, 2, 0], [1, 3, 0], [1, 4, 0], [2, 0, 0], [2, 1, 0], [2, 2, 0], [2, 3, 0], 

  [2, 4, 0], [3, 0, 0], [3, 1, 0], [3, 2, 0], [3, 3, 0], [3, 4, 0], [4, 0, 0], 

  [4, 1, 0], [4, 2, 0], [4, 3, 0], [4, 4, 0], [0, 0, 2]]
> eQ1_surfp:=[op(eQ1),1];
[1, 1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 

  1, 1]
> TQ1_surf:={}:
> for simp_nr from 1 to nops(TQ1) do
>    TQ1_surf:={op(TQ1_surf), [op(TQ1[simp_nr]), 26]}
> end: TQ1_surf;
{[1, 2, 10, 26], [1, 6, 7, 26], [1, 7, 20, 26], [1, 8, 9, 26], [1, 8, 14, 26], 

  [1, 9, 10, 26], [1, 14, 20, 26], [2, 3, 10, 26], [3, 4, 10, 26], 

  [4, 5, 10, 26], [6, 7, 13, 26], [6, 11, 12, 26], [6, 12, 25, 26], 

  [6, 13, 19, 26], [6, 19, 25, 26], [7, 13, 20, 26], [8, 9, 15, 26], 

  [8, 14, 15, 26], [9, 10, 15, 26], [11, 12, 18, 26], [11, 16, 17, 26], 

  [11, 17, 18, 26], [12, 18, 25, 26], [13, 19, 20, 26], [14, 15, 20, 26], 

  [16, 17, 25, 26], [16, 21, 22, 26], [16, 22, 23, 26], [16, 23, 24, 26], 

  [16, 24, 25, 26], [17, 18, 25, 26], [19, 20, 25, 26]}
> hQ1_surfp:=HS_Homology(Q_surf, TQ1_surf, eQ1_surfp, check_T = false, Euler = true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
         [chi = 16, h0 = _OBJ(MODZ, 9, []), h1 = _OBJ(MODZ, 2, []), 

           h2 = _OBJ(MODZ, 9, [])]
> hhQ1_surfp:=HS_NrOfComp(Q_surf, TQ1_surf, eQ1_surfp, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      9
> 
> eQ1_surfn:=[op(eQ1),-1];
[1, 1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 

  1, -1]
> hQ1_surfn:=HS_Homology(Q_surf, TQ1_surf, eQ1_surfn, check_T = false, Euler = true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
        [chi = -16, h0 = _OBJ(MODZ, 1, []), h1 = _OBJ(MODZ, 18, []), 

          h2 = _OBJ(MODZ, 1, [])]
> hhQ1_surfn:=HS_NrOfComp(Q_surf, TQ1_surf, eQ1_surfn, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      1
> 
> P_surf:=[]:
> for pkt_nr from 1 to nops(P) do
>    P_surf:=[op(P_surf), [op(P[pkt_nr]),0]]
> end do:
> P_surf:=[op(P_surf), [0,0,2]];
[[0, 0, 0], [0, 1, 0], [0, 2, 0], [0, 3, 0], [0, 4, 0], [0, 5, 0], [0, 6, 0], 

  [1, 0, 0], [1, 1, 0], [1, 2, 0], [1, 3, 0], [1, 4, 0], [1, 5, 0], [2, 0, 0], 

  [2, 1, 0], [2, 2, 0], [2, 3, 0], [2, 4, 0], [3, 0, 0], [3, 1, 0], [3, 2, 0], 

  [3, 3, 0], [4, 0, 0], [4, 1, 0], [4, 2, 0], [5, 0, 0], [5, 1, 0], [6, 0, 0], 

  [0, 0, 2]]
> TP1_surf:={}:
> for simp_nr from 1 to nops(TP1) do
>    TP1_surf:={op(TP1_surf), [op(TP1[simp_nr]), 29]}
> end: TP1_surf;
{[1, 2, 8, 29], [2, 3, 9, 29], [2, 8, 9, 29], [3, 4, 10, 29], [3, 9, 10, 29], 

  [4, 5, 10, 29], [5, 6, 11, 29], [5, 10, 15, 29], [5, 11, 15, 29], 

  [6, 7, 12, 29], [6, 11, 19, 29], [6, 12, 16, 29], [6, 16, 19, 29], 

  [7, 12, 13, 29], [8, 9, 14, 29], [9, 10, 14, 29], [10, 14, 15, 29], 

  [11, 15, 19, 29], [12, 13, 17, 29], [12, 16, 17, 29], [13, 17, 18, 29], 

  [14, 15, 19, 29], [16, 17, 20, 29], [16, 19, 20, 29], [17, 18, 21, 29], 

  [17, 20, 21, 29], [18, 21, 22, 29], [19, 20, 23, 29], [20, 21, 23, 29], 

  [21, 22, 25, 29], [21, 23, 24, 29], [21, 24, 25, 29], [23, 24, 26, 29], 

  [24, 25, 26, 29], [25, 26, 27, 29], [26, 27, 28, 29]}
> eP1_surfp:=[op(eP1),1];
[1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 

  -1, -1, 1, 1, -1, 1]
> hP1_surfp:=HS_Homology(P_surf, TP1_surf, eP1_surfp, check_T = false, Euler = true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
         [chi = 0, h0 = _OBJ(MODZ, 6, []), h1 = _OBJ(MODZ, 12, []), 

           h2 = _OBJ(MODZ, 6, [])]
> hhP1_surfp:=HS_NrOfComp(P_surf, TP1_surf, eP1_surfp, check_T = false, Euler = true);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                              [chi = 0, Nr = 6]
> eP1_surfn:=[op(eP1),-1];
[1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 

  -1, -1, 1, 1, -1, -1]
> hP1_surfn:=HS_Homology(P_surf, TP1_surf, eP1_surfn);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
  [h0 = _OBJ(MODZ, 6, []), h1 = _OBJ(MODZ, 10, []), h2 = _OBJ(MODZ, 6, [])]
> hhP1_surfn:=HS_NrOfComp(P_surf, TP1_surf, eP1_surfn, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      6
> 
> P_surf_var:=[op(1..-2,P_surf),[2,2,2]]:
> hP1_surfp_var:=HS_Homology(P_surf_var, TP1_surf, eP1_surfp, dim_range=1..1, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
                          [h1 = _OBJ(MODZ, 12, [])]
> hhP1_surfp_var:=HS_NrOfComp(P_surf_var, TP1_surf, eP1_surfp, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                             hhP1_surfp_var := 6
> hP1_surfn_var:=HS_Homology(P_surf_var, TP1_surf, eP1_surfn, dim_range=0..0, coeff=integer);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
                          [h0 = _OBJ(MODZ, 6, [])]
> hhP1_surfn_var:=HS_NrOfComp(P_surf_var, TP1_surf, eP1_surfn, Euler = true);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                              [chi = 2, Nr = 6]
> Q_surf_var:=[op(1..-2,Q_surf),[2,2,2]]:
> hQ1_surfp_var:=HS_Homology(Q_surf_var, TQ1_surf, eQ1_surfp);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
  [h0 = _OBJ(MODZ, 9, []), h1 = _OBJ(MODZ, 2, []), h2 = _OBJ(MODZ, 9, [])]
> hhQ1_surfp_var:=HS_NrOfComp(Q_surf_var, TQ1_surf, eQ1_surfp, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      9
> hQ1_surfn_var:=HS_Homology(Q_surf_var, TQ1_surf, eQ1_surfn);
                        "Checking the triangulation"
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
  [h0 = _OBJ(MODZ, 1, []), h1 = _OBJ(MODZ, 18, []), h2 = _OBJ(MODZ, 1, [])]
> hhQ1_surfn_var:=HS_NrOfComp(Q_surf_var, TQ1_surf, eQ1_surfn, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      1
> 
> W3:=[[0,0,0],[-1,0,0],[1,0,0],[0,-1,0],[0,1,0],[0,0,-1],[0,0,1],[-1,-1,0],[-1,1,0],[1,-1,0],[1,1,0],[-1,0,-1],[-1,0,1],[1,0,-1],[1,0,1],[0,-1,-1],[0,-1,1],[0,1,-1],[0,1,1],[-1,-1,-1],[-1,-1,1],[-1,1,-1],[-1,1,1],[1,-1,-1],[1,-1,1],[1,1,-1], [1,1,1]];
> 
[[0, 0, 0], [-1, 0, 0], [1, 0, 0], [0, -1, 0], [0, 1, 0], [0, 0, -1], 

  [0, 0, 1], [-1, -1, 0], [-1, 1, 0], [1, -1, 0], [1, 1, 0], [-1, 0, -1], 

  [-1, 0, 1], [1, 0, -1], [1, 0, 1], [0, -1, -1], [0, -1, 1], [0, 1, -1], 

  [0, 1, 1], [-1, -1, -1], [-1, -1, 1], [-1, 1, -1], [-1, 1, 1], [1, -1, -1], 

  [1, -1, 1], [1, 1, -1], [1, 1, 1]]
> nops(W3);
                                     27
> TW3a:={};
> for x1 from -3 to 3 do
>    for x2 from -2 to 2 do
>       for x3 from -1 to 1 do
>          v0:=[0,0,0];
>          Simplex:={v0}:
>          v1:=v0:
>          if x1 < 0 then
>             v1[-x1]:=-1
>          elif x1 > 0 then
>             v1[x1]:=1
>          else next
>          end if:
>          Simplex:=Simplex union {v1}:
>          v2:=v1:
>          if x2 < 0 then
>             k:=0:
>             for j from 1 to 3 while k < -x2 do
>                if v1[j] = 0 then k:=k+1 end if:
>             end do:
>             v2[j-1]:=-1
>          elif x2 > 0 then
>            k:=0:
>            for j from 1 to 3 while k < x2 do
>               if v1[j] = 0 then k:=k+1 end if:
>            end do:
>            v2[j-1]:=1
>          else next
>          end if:
>          Simplex:=Simplex union {v2}:
>          v3:=v2:
>          if x3 < 0 then
>             k:=0:
>             for j from 1 to 3 while k < -x3 do
>                if v2[j] = 0 then k:=k+1 end if:
>             end do:
>             v3[j-1]:=-1
>          elif x3 > 0 then
>            k:=0:
>            for j from 1 to 3 while k < x3 do
>               if v2[j] = 0 then k:=k+1 end if:
>            end do:
>            v3[j-1]:=1
>          else next
>          end if:
>          Simplex:=Simplex union {v3}:
>          TW3a:=TW3a union {Simplex}:
>       end do:
>    end do:
> end do:
                                     {}
> TW3a;
           {{[-1, -1, -1], [-1, -1, 0], [-1, 0, 0], [0, 0, 0]}, 

             {[-1, -1, -1], [-1, -1, 0], [0, -1, 0], [0, 0, 0]}, 

             {[-1, -1, -1], [-1, 0, -1], [-1, 0, 0], [0, 0, 0]}, 

             {[-1, -1, -1], [-1, 0, -1], [0, 0, -1], [0, 0, 0]}, 

             {[-1, -1, -1], [0, -1, -1], [0, -1, 0], [0, 0, 0]}, 

             {[-1, -1, -1], [0, -1, -1], [0, 0, -1], [0, 0, 0]}, 

             {[-1, -1, 0], [-1, -1, 1], [-1, 0, 0], [0, 0, 0]}, 

             {[-1, -1, 0], [-1, -1, 1], [0, -1, 0], [0, 0, 0]}, 

             {[-1, -1, 1], [-1, 0, 0], [-1, 0, 1], [0, 0, 0]}, 

             {[-1, -1, 1], [-1, 0, 1], [0, 0, 0], [0, 0, 1]}, 

             {[-1, -1, 1], [0, -1, 0], [0, -1, 1], [0, 0, 0]}, 

             {[-1, -1, 1], [0, -1, 1], [0, 0, 0], [0, 0, 1]}, 

             {[-1, 0, -1], [-1, 0, 0], [-1, 1, -1], [0, 0, 0]}, 

             {[-1, 0, -1], [-1, 1, -1], [0, 0, -1], [0, 0, 0]}, 

             {[-1, 0, 0], [-1, 0, 1], [-1, 1, 1], [0, 0, 0]}, 

             {[-1, 0, 0], [-1, 1, -1], [-1, 1, 0], [0, 0, 0]}, 

             {[-1, 0, 0], [-1, 1, 0], [-1, 1, 1], [0, 0, 0]}, 

             {[-1, 0, 1], [-1, 1, 1], [0, 0, 0], [0, 0, 1]}, 

             {[-1, 1, -1], [-1, 1, 0], [0, 0, 0], [0, 1, 0]}, 

             {[-1, 1, -1], [0, 0, -1], [0, 0, 0], [0, 1, -1]}, 

             {[-1, 1, -1], [0, 0, 0], [0, 1, -1], [0, 1, 0]}, 

             {[-1, 1, 0], [-1, 1, 1], [0, 0, 0], [0, 1, 0]}, 

             {[-1, 1, 1], [0, 0, 0], [0, 0, 1], [0, 1, 1]}, 

             {[-1, 1, 1], [0, 0, 0], [0, 1, 0], [0, 1, 1]}, 

             {[0, -1, -1], [0, -1, 0], [0, 0, 0], [1, -1, -1]}, 

             {[0, -1, -1], [0, 0, -1], [0, 0, 0], [1, -1, -1]}, 

             {[0, -1, 0], [0, -1, 1], [0, 0, 0], [1, -1, 1]}, 

             {[0, -1, 0], [0, 0, 0], [1, -1, -1], [1, -1, 0]}, 

             {[0, -1, 0], [0, 0, 0], [1, -1, 0], [1, -1, 1]}, 

             {[0, -1, 1], [0, 0, 0], [0, 0, 1], [1, -1, 1]}, 

             {[0, 0, -1], [0, 0, 0], [0, 1, -1], [1, 1, -1]}, 

             {[0, 0, -1], [0, 0, 0], [1, -1, -1], [1, 0, -1]}, 

             {[0, 0, -1], [0, 0, 0], [1, 0, -1], [1, 1, -1]}, 

             {[0, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1]}, 

             {[0, 0, 0], [0, 0, 1], [1, -1, 1], [1, 0, 1]}, 

             {[0, 0, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1]}, 

             {[0, 0, 0], [0, 1, -1], [0, 1, 0], [1, 1, -1]}, 

             {[0, 0, 0], [0, 1, 0], [0, 1, 1], [1, 1, 1]}, 

             {[0, 0, 0], [0, 1, 0], [1, 1, -1], [1, 1, 0]}, 

             {[0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]}, 

             {[0, 0, 0], [1, -1, -1], [1, -1, 0], [1, 0, 0]}, 

             {[0, 0, 0], [1, -1, -1], [1, 0, -1], [1, 0, 0]}, 

             {[0, 0, 0], [1, -1, 0], [1, -1, 1], [1, 0, 0]}, 

             {[0, 0, 0], [1, -1, 1], [1, 0, 0], [1, 0, 1]}, 

             {[0, 0, 0], [1, 0, -1], [1, 0, 0], [1, 1, -1]}, 

             {[0, 0, 0], [1, 0, 0], [1, 0, 1], [1, 1, 1]}, 

             {[0, 0, 0], [1, 0, 0], [1, 1, -1], [1, 1, 0]}, 

             {[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]}}
> nops(TW3a);
                                     48
> TW3:={}:
> for Simplex in TW3a do
>    Smp:={};
>    for vertex in Simplex do
>       for j from 1 to nops(W3) while vertex <> W3[j] do end do:
>       Smp:=Smp union {j}:
>    end do:
>    TW3:=TW3 union {Smp}:
> end do:
> TW3;
{{1, 2, 8, 20}, {1, 2, 8, 21}, {1, 2, 9, 22}, {1, 2, 9, 23}, {1, 2, 12, 20}, 

  {1, 2, 12, 22}, {1, 2, 13, 21}, {1, 2, 13, 23}, {1, 3, 10, 24}, 

  {1, 3, 10, 25}, {1, 3, 11, 26}, {1, 3, 11, 27}, {1, 3, 14, 24}, 

  {1, 3, 14, 26}, {1, 3, 15, 25}, {1, 3, 15, 27}, {1, 4, 8, 20}, 

  {1, 4, 8, 21}, {1, 4, 10, 24}, {1, 4, 10, 25}, {1, 4, 16, 20}, 

  {1, 4, 16, 24}, {1, 4, 17, 21}, {1, 4, 17, 25}, {1, 5, 9, 22}, 

  {1, 5, 9, 23}, {1, 5, 11, 26}, {1, 5, 11, 27}, {1, 5, 18, 22}, 

  {1, 5, 18, 26}, {1, 5, 19, 23}, {1, 5, 19, 27}, {1, 6, 12, 20}, 

  {1, 6, 12, 22}, {1, 6, 14, 24}, {1, 6, 14, 26}, {1, 6, 16, 20}, 

  {1, 6, 16, 24}, {1, 6, 18, 22}, {1, 6, 18, 26}, {1, 7, 13, 21}, 

  {1, 7, 13, 23}, {1, 7, 15, 25}, {1, 7, 15, 27}, {1, 7, 17, 21}, 

  {1, 7, 17, 25}, {1, 7, 19, 23}, {1, 7, 19, 27}}
> checkTriangulation(W3, TW3);
                                    true
> eW3:=[seq(1, i=1..27)];
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 

  1]
> eW3[1]:=-1;
                                     -1
> hW3:=HS_Homology(W3, TW3, [seq(1, i=1..27)], Euler=true, check_T=false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
        [chi = -16, h0 = _OBJ(MODZ, 1, []), h1 = _OBJ(MODZ, 18, []), 

          h2 = _OBJ(MODZ, 1, [])]
> nrChW3:=HS_NrOfComp(W3, TW3, eW3, Euler=true, check_T=false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                             [chi = -16, Nr = 2]
> 
> W0:=[[]]:
> for n from 1 to 4 do
>    W||n:=[]:
>    for pt in W||(n-1) do
>       W||n:=[op(W||n), [op(pt), 1],
>                        [op(pt), 0],
>                        [op(pt), -1] ]
>    end do:
> end do: 
> W4;
[[1, 1, 1, 1], [1, 1, 1, 0], [1, 1, 1, -1], [1, 1, 0, 1], [1, 1, 0, 0], 

  [1, 1, 0, -1], [1, 1, -1, 1], [1, 1, -1, 0], [1, 1, -1, -1], [1, 0, 1, 1], 

  [1, 0, 1, 0], [1, 0, 1, -1], [1, 0, 0, 1], [1, 0, 0, 0], [1, 0, 0, -1], 

  [1, 0, -1, 1], [1, 0, -1, 0], [1, 0, -1, -1], [1, -1, 1, 1], [1, -1, 1, 0], 

  [1, -1, 1, -1], [1, -1, 0, 1], [1, -1, 0, 0], [1, -1, 0, -1], 

  [1, -1, -1, 1], [1, -1, -1, 0], [1, -1, -1, -1], [0, 1, 1, 1], [0, 1, 1, 0], 

  [0, 1, 1, -1], [0, 1, 0, 1], [0, 1, 0, 0], [0, 1, 0, -1], [0, 1, -1, 1], 

  [0, 1, -1, 0], [0, 1, -1, -1], [0, 0, 1, 1], [0, 0, 1, 0], [0, 0, 1, -1], 

  [0, 0, 0, 1], [0, 0, 0, 0], [0, 0, 0, -1], [0, 0, -1, 1], [0, 0, -1, 0], 

  [0, 0, -1, -1], [0, -1, 1, 1], [0, -1, 1, 0], [0, -1, 1, -1], [0, -1, 0, 1], 

  [0, -1, 0, 0], [0, -1, 0, -1], [0, -1, -1, 1], [0, -1, -1, 0], 

  [0, -1, -1, -1], [-1, 1, 1, 1], [-1, 1, 1, 0], [-1, 1, 1, -1], 

  [-1, 1, 0, 1], [-1, 1, 0, 0], [-1, 1, 0, -1], [-1, 1, -1, 1], 

  [-1, 1, -1, 0], [-1, 1, -1, -1], [-1, 0, 1, 1], [-1, 0, 1, 0], 

  [-1, 0, 1, -1], [-1, 0, 0, 1], [-1, 0, 0, 0], [-1, 0, 0, -1], 

  [-1, 0, -1, 1], [-1, 0, -1, 0], [-1, 0, -1, -1], [-1, -1, 1, 1], 

  [-1, -1, 1, 0], [-1, -1, 1, -1], [-1, -1, 0, 1], [-1, -1, 0, 0], 

  [-1, -1, 0, -1], [-1, -1, -1, 1], [-1, -1, -1, 0], [-1, -1, -1, -1]]
> nops(W4);
                                     81
>  Verts:={[]}:
>    for i from 1 to 4 do
>       Verts:={seq([1,op(k)], k in Verts), seq([-1,op(k)], k in Verts)}
>    end do:
> Verts;
    {[-1, -1, -1, -1], [-1, -1, -1, 1], [-1, -1, 1, -1], [-1, -1, 1, 1], 

      [-1, 1, -1, -1], [-1, 1, -1, 1], [-1, 1, 1, -1], [-1, 1, 1, 1], 

      [1, -1, -1, -1], [1, -1, -1, 1], [1, -1, 1, -1], [1, -1, 1, 1], 

      [1, 1, -1, -1], [1, 1, -1, 1], [1, 1, 1, -1], [1, 1, 1, 1]}
> l_TW0:={seq([Verts[i]], i=1..nops(Verts))}:
> for n from 1 to 4 do
>    l_TW||n:={}:
>    for simp in l_TW||(n-1) do
>       v:=simp[-1]:
>       for k from 1 to nops(v) do
>          if v[k] <> 0 then
>             simp_new:=[op(simp), subsop(k=0, v)]:
>             l_TW||n:=l_TW||n union {simp_new}
>          end if
>       end do
>    end do
> end do:
> nops(l_TW4);
                                     384
> TW4:={}:
> for Simplex in l_TW4 do
>    Smp:={};
>    for vertex in Simplex do
>       for j from 1 to nops(W4) while vertex <> W4[j] do end do:
>       Smp:=Smp union {j}:
>    end do:
>    TW4:=TW4 union {Smp}:
> end do:
> TW4[1];
                              {1, 2, 5, 14, 41}
> checkTriangulation(W4, TW4);
                                    true
> eW4:=[seq(1, i=1..nops(W4))]:
> hW4:=HS_Homology(W4, TW4, eW4, check_T = false, Euler = false, coeff=mod2, dim_range=1..1);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                           "Calculating the ranks"
> nrCW4:=HS_NrOfComp(W4, TW4, eW4, check_T = false, Euler = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                      1
> 
> 
> eW4_mult:=proc(i::posint)::procedure;
> 
> local d, k, l, vorz_struct, vorz, LattHom2_store;
> 
> global W4, LattHom2;
> 
> d:=4;
> if nops(LattHom2[1]) <> d then error "Falsches LattHom2" end if;
> if i > 2^(d+1) then return [] end if;
> LattHom2_store:=LattHom2;
> LattHom2:=[seq([1,op(k)], k in LattHom2), seq([-1,op(k)], k in LattHom2)];
> vorz_struct:=Nr2Hom(i);
> vorz:=[seq(0, l=1..3^d)];
> for l from 1 to 3^d do
>    vorz[l]:=vorz_struct[nops(select(x -> x=0, W4[l]))+1]
> end do;
> LattHom2:=LattHom2_store;
> vorz
> end proc:
> nrCW4_mult:=HS_NrOfComp(W4, TW4, eW4_mult, check_T = false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                           "Calculating the rest"
                                   {1, 2}
> eW4_mult(32);
[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 

  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 

  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 

  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 

  -1, -1, -1, -1, -1]
> stopat(eW4_mult,1);
                                 [eW4_mult]
> LattHom2:=[[]];
>    for i from 1 to 4 do
>       LattHom2:=[seq([1,op(k)], k in LattHom2), seq([-1,op(k)], k in LattHom2)]
>    end do:
                                    [[]]
> nops(LattHom2[1]);
                                      5
> C3:=[[0,0,0]]:
> for i from 1 to 3 do
>    for j in [1,-1] do
>       pt:=[seq(0, k=1..3)]:
>       pt[i]:=j:
>       C3:=[op(C3), pt]
>    end do
> end do: 
> C3;
[[0, 0, 0], [1, 0, 0], [-1, 0, 0], [0, 1, 0], [0, -1, 0], [0, 0, 1], [0, 0, -1]

  ]
> TC3:={}:
> for simp_nr from 1 to 8 do
>    simp:={1}:
>    n:=simp_nr:
>    for k from 3 to 1 by -1 do
>       nmod:=irem(n,2, 'n'):
>       v:=2*k+nmod:
>       simp:=simp union {v}
>    end do:
> TC3:=TC3 union {simp}
> end do:
> TC3;
   {{1, 2, 4, 6}, {1, 2, 4, 7}, {1, 2, 5, 6}, {1, 2, 5, 7}, {1, 3, 4, 6}, 

     {1, 3, 4, 7}, {1, 3, 5, 6}, {1, 3, 5, 7}}
> checkTriangulation(C3, TC3);
                                    true
> eC3:=[seq(1, i=1..7)];
                            [1, 1, 1, 1, 1, 1, 1]
> hC3:=HS_Homology(C3, TC3, eC3, Euler=true, check_T =false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
        [chi = -4, h0 = _OBJ(MODZ, 2, []), h1 = _OBJ(MODZ, 10, []), 

          h2 = _OBJ(MODZ, 4, [])]
> 
> C4:=[[0,0,0,0]]:
> for i from 1 to 4 do
>    for j in [1,-1] do
>       pt:=[seq(0, k=1..4)]:
>       pt[i]:=j:
>       C4:=[op(C4), pt]
>    end do
> end do: 
> C4;
  [[0, 0, 0, 0], [1, 0, 0, 0], [-1, 0, 0, 0], [0, 1, 0, 0], [0, -1, 0, 0], 

    [0, 0, 1, 0], [0, 0, -1, 0], [0, 0, 0, 1], [0, 0, 0, -1]]
> TC4:={}:
> for simp_nr from 1 to 16 do
>    simp:={1}:
>    n:=simp_nr:
>    for k from 4 to 1 by -1 do
>       nmod:=irem(n,2, 'n'):
>       v:=2*k+nmod:
>       simp:=simp union {v}
>    end do:
> TC4:=TC4 union {simp}
> end do:
> TC4;
   {{1, 2, 4, 6, 8}, {1, 2, 4, 6, 9}, {1, 2, 4, 7, 8}, {1, 2, 4, 7, 9}, 

     {1, 2, 5, 6, 8}, {1, 2, 5, 6, 9}, {1, 2, 5, 7, 8}, {1, 2, 5, 7, 9}, 

     {1, 3, 4, 6, 8}, {1, 3, 4, 6, 9}, {1, 3, 4, 7, 8}, {1, 3, 4, 7, 9}, 

     {1, 3, 5, 6, 8}, {1, 3, 5, 6, 9}, {1, 3, 5, 7, 8}, {1, 3, 5, 7, 9}}
> checkTriangulation(C4, TC4);
                                    true
> eC4:=[seq(1, i=1..9)];
                         [1, 1, 1, 1, 1, 1, 1, 1, 1]
> eC4[2]:=-1;eC4[4]:=-1;eC4[6]:=-1;eC4[8]:=-1;
                                     -1
                                     -1
                                     -1
                                     -1
> 
> hC4:=HS_Homology(C4, TC4, eC4, Euler=true, check_T =false, coeff=integer);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
        [chi = 24, h0 = _OBJ(MODZ, 2, []), h1 = _OBJ(MODZ, 4, [2]), 

          h2 = _OBJ(MODZ, 34, []), h3 = _OBJ(MODZ, 8, [])]
> nrC_C4:=HS_NrOfComp(C4, TC4, eC4, check_T =false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                                 nrC_C4 := 2
> 
> 
> Tet4:=[[0,0,0,0], seq(subsop(i=1, [0,0,0,0]), i=1..4), [-1,-1,-1,-1]];
   [[0, 0, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], 

     [-1, -1, -1, -1]]
> TTet4:={seq({1,op(subsop(i=NULL, [2,3,4,5,6]))},i=1..5)};
    {{1, 2, 3, 4, 5}, {1, 2, 3, 4, 6}, {1, 2, 3, 5, 6}, {1, 2, 4, 5, 6}, 

      {1, 3, 4, 5, 6}}
> eTet4:=[seq(1, i=1..nops(Tet4))];
                             [1, 1, 1, 1, 1, 1]
> eTet4[6]:=-1;
                                     -1
> checkTriangulation(Tet4, TTet4);
                                    true
> hTet4:=HS_Homology(Tet4, TTet4, eTet4, Euler=true, check_T =false, coeff=integer);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
       "Determine the cells of the hypersurface and their orientation"
                       "Calculating the boundary maps"
                    "Calculating the Smith normal forms"
         [chi = 0, h0 = _OBJ(MODZ, 2, []), h1 = _OBJ(MODZ, 0, [2]), 

           h2 = _OBJ(MODZ, 0, []), h3 = _OBJ(MODZ, 2, [])]
> nrCTet4:=HS_NrOfComp(Tet4, TTet4, eTet4, Euler=true, check_T=false);
              "Calculating the incidences of the triangulation"
                     "Calculating the glueing structure"
                  "Determine the cells of the hypersurface"
                   "Calculating the number of components"
                              [chi = 0, Nr = 2]
> stopat(HS_Homology, 62);
                                [HS_Homology]
> unstopat();
                                     []
>