\\ ---------------  GP code  ------------------------------------------------------------
\\ Library for the program Bianchi.gp
\\
\\ Description: Compute the quotient of Hyperbolic Space by PSL_2 of imaginary
\\ quadratic number fields 
\\
\\
\\ Author: Alexander D. Rahm
\\                  
\\ Copyright (C) 2010 by Alexander D. Rahm. 
\\ Bianchi.gp is a free software covered by the GNU General Public License.
\\ Version 2.0.4 of June 29, 2011.     
\\---------------------------------------------------------------------------------------


ComputeOrbitSpaceNot3mod4() =
{
  local( orbitSum);

  if( m > 1,
        subdivideVerticalLoop();
        prepareDivideOutAssociatedMatrices(); 
        getVertexStabilizers();
        getVertexIdentifications();
        initialize2cells();
        getVertexOrbits();
        recordCorners();
        spuriousCellFissionNot3mod4();
        divideOutAssociatedMatrices();
        sortByOrbitNumbers();
        orbitSum = recordOrbitSum();
        identify2cells( orbitSum);
        getEdges();
        computeEdgeStabilizers();
        getEdgeIdentifications();
        getEdgeOrbits();
        getOriginAndEndNumbers();
        computeEdgeStabilizers();
  );
};
/* end of procedure ComputeOrbitSpaceNot3mod4 */


ComputeOrbitSpace3mod4() =
{
  local( orbitSum);
  if( m > 3, /* Get the quotient by the group action. */

        subdivideVerticalLoop();
        getVertexStabilizers();                 
        getVertexIdentifications();
        getVertexOrbits();
        initialize2cells();
        recordCorners();
        divideOutAssociatedMatrices();
        spuriousCellFission();
        sortByOrbitNumbers(); 
        orbitSum = recordOrbitSum();
        getEdges(); 
        computeEdgeStabilizers();
        getEdgeIdentifications();
        getEdgeOrbits();
        identify2cells( orbitSum);
        getOriginAndEndNumbers();
        computeEdgeStabilizers();
  );
};
/* end of procedure ComputeOrbitSpace3mod4 */



getEdges() =
{
  local( edge, auxiliaryedgesOf2cellList:list, existingEdgeNumber, overchargeCounter);   
  local( cornersOf2cell_j );
  overchargeCounter = 0;
  edgesList = listcreate(2*numberOfLines);
  listkill(edgesList:list);
  /* Intersect the sets pointsOfLine[ twoCellSupport] and cornersOf2cell for every 2-cell. */
  /* Write the edges in this intersection, consisting of their origin and end point, */
  /* into edgesList. */
  for (j = 1, numberOf2cells, 

      /* Create the set of coordinates of elements of cornersOf2cell[j]. */
      cornersOf2cell_j = vector( length(cornersOf2cell[j]));
      for( s = 1, length( cornersOf2cell[j]),

                cornersOf2cell_j[s] = eval(totalPointSet[ cornersOf2cell[j][s]]); 
      );
      for( s = 1, length( pointsOfLine[ twoCellSupport:list[j]]),

         edge = setintersect( Set( pointsOfLine[ twoCellSupport:list[j]][s]),
                                Set( cornersOf2cell_j));

         if(length(edge) == 2, 

                listput(edgesList:list,edge);   
         );
         if(length(edge) > 2,
                overchargeCounter++;
         );
      ); 
  );
  print("There are ",overchargeCounter," edges which are overcharged with vertices.");

  /* Kill double entries with the parameter "1" in the following command: */
  listsort(edgesList:list, 1);

  /* In order to find the entries of edgesList with setsearch, */
  /* give edgesList the ordering of Set(edgesList). */ 
  for ( k = 1, length(edgesList:list),

        /* Replace the k'th entry of edgesList with the k'th entry of Set(edgesList). */
        listput( edgesList:list, eval(Set(edgesList:list)[k]), k);
  );

  edgesOf2cell = vector(numberOf2cells);
  /* Record into the lists edgesOf2cell[j] the numbers in the sorted list edgesList */
  /* of edges lying in the hemisphere number j. */
  for (j = 1, numberOf2cells, 

      auxiliaryedgesOf2cellList = listcreate( length( linesOfSphere[ twoCellSupport:list[j]]) +2);

      /* Create the set of coordinates of elements of cornersOf2cell[j]. */
      cornersOf2cell_j = vector( length(cornersOf2cell[j]));
      for( s = 1, length( cornersOf2cell[j]),

                cornersOf2cell_j[s] = eval(totalPointSet[ cornersOf2cell[j][s]]); 
      );

      for( s = 1, length( pointsOfLine[ twoCellSupport:list[j]]),

         edge = setintersect( Set( pointsOfLine[ twoCellSupport:list[j]][s]), Set( cornersOf2cell_j));

         if(length(edge) == 2, 

            /* As we have given edgesList the ordering of Set(edgesList), */
            /* we may use the number */
            existingEdgeNumber = setsearch( Set(edgesList:list), edge);
            if(existingEdgeNumber,
                
                listput(auxiliaryedgesOf2cellList, existingEdgeNumber);
            );  
         );
      );
      edgesOf2cell[j] = auxiliaryedgesOf2cellList;
      listkill( auxiliaryedgesOf2cellList);
  );
};
{addhelp(getEdges, "Computes the edges for existing global data on the lines and corners.");}
/* end of procedure getEdges. */



computeEdgeStabilizers() =
{
  /* Compute the stabilizers of all edges in edgesList */
        /* and record them as the global vector edgeStabilizer. */
  local(OriginStabilizer,EndStabilizer, z, rsquare);
  local( edgeOrigin, edgeEnd);

  edgeStabilizer = vector( length( edgesList:list));

  for ( j = 1, length( edgesList:list),

        edgeOrigin = eval(edgesList:list[j][1]);
        edgeEnd = eval(edgesList:list[j][2]);

        /* edgeStabilizer[j] 
                = intersection( stabilizer(edgeOrigin), stabilizer(edgeEnd)); */

        /* get OriginStabilizer */
        z = component(edgeOrigin,1);
        rsquare = component(edgeOrigin,2);

        if ( rsquare == 0,  /* vertex is a singular point */
        
          OriginStabilizer = Set ( computeSingularStabilizer(z));

        , /* else vertex is in the interior of Hyperbolic Space */
          if ( mIs3mod4,

                  OriginStabilizer = Set ( computeVertexStabilizer3mod4(z, rsquare));

          , /* else the one or two modulo 4 function call */

                  OriginStabilizer = Set ( computeVertexStabilizer(z, rsquare));
          );
        );

        /* get EndStabilizer */
        z = component(edgeEnd,1);
        rsquare = component(edgeEnd,2);

        if ( rsquare == 0,  /* vertex is a singular point */
        
          EndStabilizer = Set ( computeSingularStabilizer(z));

        , /* else vertex is in the interior of Hyperbolic Space */
          if ( mIs3mod4,

                  EndStabilizer = Set ( computeVertexStabilizer3mod4(z, rsquare));

          , /* else the one or two modulo 4 function call */

                  EndStabilizer = Set ( computeVertexStabilizer(z, rsquare));
          );
        );

        edgeStabilizer[j] = setintersect(OriginStabilizer,EndStabilizer);
  );
};
/* end of procedure computeEdgeStabilizers */




getVertexStabilizers()=
{
  local( z, rsquare);
  
  if (m == 3, 
        print("Stabilizers not completely calculated because of additional units in the ", 
                        "ring of integers of Q(sqrt(-3)).")
  );
  stabilizer = vector(length(totalPointSet));

  for( j = 1, length( totalPointSet),

      z = component(eval(totalPointSet[j]),1);
      rsquare = component(eval(totalPointSet[j]),2);

      if ( rsquare == 0,  /* vertex is a singular point */
        
          stabilizer[j] = computeSingularStabilizer(z);

        , /* else vertex is in the interior of Hyperbolic Space */
          if ( mIs3mod4,

                  stabilizer[j] = computeVertexStabilizer3mod4( z, rsquare);

          , /* else the one or two modulo 4 function call */

                  stabilizer[j] = computeVertexStabilizer( z, rsquare);
          );
      );       
  );
  print(length(totalPointSet)," vertex stabilizers computed.");
};
/* end of procedure getVertexStabilizers. */




getVertexIdentifications() = 
{
  /* The set identificationMatrices[index] shall contain the matrices transporting */
    /* one point into another in totalPointSet. */ 
  
  local( vectorOfIdentifications, singular_r: list, singular_s: list, point_r, point_s, vertexIdentification, auxiliaryList);

  singular_r = listcreate( 4*length(singularpoint)^2);
  singular_s = listcreate( 4*length(singularpoint)^2);

  identificationMatrices = listcreate( 2*length(totalPointSet)*maxNeighboursNumber);
  transportFrom = vector( length(totalPointSet) );
  equivalentVertices = vector( length(totalPointSet) );
  for( r = 1, length(totalPointSet),
    transportFrom[r] = listcreate(2*maxNeighboursNumber);
    equivalentVertices[r] = listcreate(2*maxNeighboursNumber);
    listput( identificationMatrices, Set(stabilizer[r]) );
    transportFrom[r] = sublistPut(transportFrom[r], length(identificationMatrices));
    equivalentVertices[r] = sublistPut( equivalentVertices[r], r);      
  );
  /* Consider each pair of vertices, one after another. */
  for( r = 1, length(totalPointSet),

    point_r = eval(totalPointSet[r]);

    for( s = r+1, length(totalPointSet),

      point_s = eval(totalPointSet[s]);

      /* Check that the stabilizers of both vertices have the same isomorphy type. */
      /*It suffices to check the cardinal, which classifies the isomorphy type here. */ 
      if( length( stabilizer[r]:list) == length( stabilizer[s]:list),

        /*Check that both point_r and point_s are in the interior of hyperbolic space.*/
        if( component( point_r, 2) > 0 && component( point_s, 2) > 0,

          /*Call the function getIdentificationMatrices that computes the 
                identification matrices between the pair.*/
          if( mIs3mod4,

             vertexIdentification = getIdentificationMatrices3mod4(
                                        point_r, point_s, length( stabilizer[r]:list));
             if( length(vertexIdentification) > 0, 
                        listput( identificationMatrices, Set(vertexIdentification) );
                transportFrom[r] = sublistPut(transportFrom[r], length(identificationMatrices));
                equivalentVertices[r] = sublistPut( equivalentVertices[r], s);     
             /* The following operation can be replaced taking the inverse matrices.*/
                        vertexIdentification =  getIdentificationMatrices3mod4(
                                        point_s, point_r, length( stabilizer[r]:list)) ;
                        listput( identificationMatrices, Set(vertexIdentification) );
                transportFrom[s] = sublistPut( transportFrom[s], length(identificationMatrices));
                equivalentVertices[s] = sublistPut( equivalentVertices[s], r);     
             );
          ,/* else m not congruent 3 modulo 4. */            
             vertexIdentification = getIdentificationMatrices(
                                        point_r, point_s, length( stabilizer[r]:list)) ;
             if( length(vertexIdentification) > 0, 
                        listput( identificationMatrices, Set(vertexIdentification) );
                transportFrom[r] = sublistPut( transportFrom[r], length(identificationMatrices));
                equivalentVertices[r] = sublistPut( equivalentVertices[r], s);     
             /* The following operation can be replaced taking the inverse matrices.*/
                        vertexIdentification =  getIdentificationMatrices(
                                        point_s, point_r, length( stabilizer[r]:list)) ;
                        listput( identificationMatrices, Set(vertexIdentification) );
                transportFrom[s] = sublistPut( transportFrom[s], length(identificationMatrices));
                equivalentVertices[s] = sublistPut( equivalentVertices[s], r);     
             );
          );

        ,/* else point_r and point_s are singular. Mark them to be treated later. */
        
          listput(  singular_r, r); listput( singular_s, s);
        );
      );
    );
  );
  vectorOfIdentifications = recordVectorOfIdentifications();


  /* Now get the identifications of singular points. */
  for( k = 1, length( singular_r),
    vertexIdentification = getSingularIdentifications(
                                singular_r[k], singular_s[k], vectorOfIdentifications);
    if( length(vertexIdentification) > 0,
        listput( identificationMatrices, Set(vertexIdentification) );
        transportFrom[singular_r[k]] = sublistPut( 
                     transportFrom[singular_r[k]], length(identificationMatrices));
        equivalentVertices[singular_r[k]] = sublistPut( 
                        equivalentVertices[singular_r[k]], singular_s[k]);      
        /* The following operation can be replaced taking the inverse matrices.*/
        vertexIdentification =  getSingularIdentifications(
                                singular_s[k], singular_r[k], vectorOfIdentifications);
        listput( identificationMatrices, Set(vertexIdentification) );
        transportFrom[singular_s[k]] = sublistPut( 
                        transportFrom[singular_s[k]], length(identificationMatrices));
        equivalentVertices[singular_s[k]] = sublistPut( 
                        equivalentVertices[singular_s[k]], singular_r[k]);      
    );
  );
  recordSingularSelfIdentifications(vectorOfIdentifications); 
  for(j = 1, 4,
    recordSingularTransitivityRelations();
  );
};
/* end of the procedure getVertexIdentifications. */


recordSingularSelfIdentifications(vectorOfIdentifications) = {
local( point_r, index);
  for( r = 1, length(totalPointSet),

    point_r = eval(totalPointSet[r]);
    if( component( point_r, 2) == 0, /* if r is singular, get its self-identifications */
        index = listSearch( equivalentVertices[r], r);
        identificationMatrices[transportFrom[r][index]] = setunion(
                Set(stabilizer[r]),
                Set(getSingularIdentifications(r, r, vectorOfIdentifications))
        );
    );
  );
};
/* end of the procedure recordSingularSelfIdentifications. */




recordSingularTransitivityRelations() = {
local( point_s, t, u, index_u_for_s);

  for( s = 1, length(totalPointSet),
    point_s = eval(totalPointSet[s]);

    if( component( point_s, 2) == 0, /* if s is singular, get its transitivity relations */

      for( index_t = 1, length( equivalentVertices[s]),

        t = equivalentVertices[s][index_t];

        if( t != s,

         for( index_u = 1, length( equivalentVertices[t]),

          u = equivalentVertices[t][index_u];

          if( u != t,

           index_u_for_s = listSearch( equivalentVertices[s], u);

           if( index_u_for_s > 0, /* if u is already recorded in the orbit of s */

            identificationMatrices[transportFrom[s][index_u_for_s]] = setunion(
                identificationMatrices[transportFrom[s][index_u_for_s]], 
                Composition(identificationMatrices[transportFrom[s][index_t]],
                            identificationMatrices[transportFrom[t][index_u]]
                )
            );
           , /* else u enters the orbit of s */
            listput( identificationMatrices, 
                        Composition(
                            identificationMatrices[transportFrom[s][index_t]],
                            identificationMatrices[transportFrom[t][index_u]]
                        )
            );
            transportFrom[s] = sublistPut( transportFrom[s], #(identificationMatrices));
            equivalentVertices[s] = sublistPut( equivalentVertices[s], u);
           );
          );
         );
        );
      );
    );
  );
};
/* end of the procedure recordSingularTransitivityRelations. */



getEdgeIdentifications() =
{
  local( formerOrigin, origin_r, end_r, origin_s, end_s, vertices_equivalent_to_end_r, vertices_equivalent_to_origin_r, o,e);

  equivalentEdges = vector( length(edgesList) );

  /* In order to identify edges, get the intersections of 
        the vertexIdentifications sets. */
  for( r = 1, length(edgesList:list),

    origin_r = setsearch(totalPointSet, component(edgesList:list[r],1));
    end_r = setsearch(totalPointSet, component(edgesList:list[r],2));
    equivalentEdges[r] = listcreate(2*maxNeighboursNumber);

    vertices_equivalent_to_end_r = equivalentVertices[end_r];
    vertices_equivalent_to_origin_r = equivalentVertices[origin_r];

    for( s = 1, length(edgesList:list),
      if( s != r,

        origin_s = setsearch(totalPointSet, component(edgesList:list[s],1));
        end_s = setsearch(totalPointSet, component(edgesList:list[s],2));
        o = listSearch( vertices_equivalent_to_origin_r, origin_s);

        if( o > 0,
          e = listSearch( vertices_equivalent_to_end_r, end_s);

          if( e > 0,

            if( length( setintersect( 
                identificationMatrices[transportFrom[origin_r][o]],
                identificationMatrices[transportFrom[ end_r][e]] )) > 0,

                    equivalentEdges[r] = sublistPut(equivalentEdges[r],s);
            );
          );
        );
        /* Now do the same again with the edge geometrically reverted: */
        formerOrigin = origin_s;
        origin_s = end_s;
        end_s = formerOrigin;
        o = listSearch( vertices_equivalent_to_origin_r, origin_s);

        if( o > 0,
          e = listSearch( vertices_equivalent_to_end_r, end_s);

          if( e > 0,

            if( length( setintersect( 
                identificationMatrices[transportFrom[origin_r][o]],
                identificationMatrices[transportFrom[ end_r][e]] )) > 0,

                    equivalentEdges[r] = sublistPut(equivalentEdges[r],s);
            );
          );
        );
      );
    );
  );
};
/* end of procedure getEdgeIdentifications */


addVertexIdentifications(j, Fraction) = 
{
  /* Write the matrix g associated to the hemisphere number k by the procedure */
  /* divideOutAssociatedMatrices() into the entries of "vertexIdentifications", */
  /* where g identifies two vertices lying both on the hemisphere number j. */
  local( gz, a,b,c,d, z, rsquare, k_index, l_index, newEdge, index);
  newEdge = vector(2);

  /* For g = [ -lambda, -(lambda^2 +1)/mu; mu, lambda] =: [a, b; c, d], */
  /* compute the images g*{vertex on hemisphere number j}. */
  a = -Lambda[j]; b = -Fraction; c = Mu[j]; d = Lambda[j];

  for( k = 1, length( pointsOfSphere[j]:list),

    /* Take the vertex (z,r) := pointsOfSphere[j][k] */
    z = component( pointsOfSphere[j]:list[k], 1);
    rsquare = component( pointsOfSphere[j]:list[k], 2);

    /* and apply g to it, using Floege's formula for [a, b; c, d]*(z,r). */
    gz =  nfeltdiv(K, nfbasistoalg(K,d) +  conj(nfbasistoalg(K, d -nfeltmul(K, c, z))),
                   c);
    if( gz != nfeltmul(K, conj(nfbasistoalg(K, d -nfeltmul(K, c, z))),
                        nfeltmul(K, a, z) -b)
           -rsquare *nfeltmul(K, conj(nfbasistoalg(K, c)), a),
                print("***Error in procedure addVertexIdentifications at the expression for Floege's operation formula."); 
    );
    /* gr = r/denominatorExpression == 1 by Floege. */
    if( 1 != norm(nfbasistoalg( K, nfeltmul(K, c, z) -d))
                + rsquare *norm(nfbasistoalg(K, c)),
        print("***Error in procedure addVertexIdentifications at the denominator expression for Floege's operation formula."); 
    );

    for( l = k, length( pointsOfSphere[j]:list),
        
        if( gz ==  component( pointsOfSphere[j]:list[l], 1),
                
                print("g* (",z," , ",rsquare," ) = ( ",
                        gz," , ",component( pointsOfSphere[j]:list[l], 2)," ).");

                /* enter g into vertexIdentifications for k,l */
                k_index = setsearch( totalPointSet, [z, rsquare]);
                l_index = setsearch( totalPointSet, [gz, rsquare]);
                print(k_index, " , ",l_index);
                index = listSearch( equivalentVertices[k_index], l_index);
                identificationMatrices[transportFrom[k_index][index]] = setunion(
                        Set( [[a,b;c,d], -[a,b;c,d]] ), 
                        identificationMatrices[transportFrom[k_index][index]]);
                /* enter g^-1 into vertexIdentifications for l,k */
                index = listSearch( equivalentVertices[l_index], k_index);
                identificationMatrices[transportFrom[l_index][index]] = setunion(
                        Set( [[d,-b;-c,a], -[d,-b;-c,a]] ), 
                        identificationMatrices[transportFrom[l_index][index]]); 
        );
    ); 
  );
};
/* end of procedure addVertexIdentifications */


insertEdge( j, newEdge)=
/* Insert newEdge on sphere number j. */
{ 
  local( RealPartOne, RealPartTwo);

  if( mIs3mod4 == 0, /* m not congruent to 3 mod 4 */

        RealPartOne = component( newEdge[1], 1);
        RealPartTwo = component( newEdge[2], 1);

  , /* else m congruent 3 modulo 4: */

        RealPartOne = component( newEdge[1], 1) -component( newEdge[1], 2)/2;
        RealPartTwo = component( newEdge[2], 1) -component( newEdge[2], 2)/2;
  );
  if( RealPartOne == RealPartTwo,

        insertVerticalEdge( j, RealPartOne);

  , /* else not a vertical edge */
        print("****Error in procedure insertEdge: The procedure insertHorizontalEdge must be reactivated and called from here");
  ); 
}; /* end of procedure insertEdge */



prepareDivideOutAssociatedMatrices() = 
{
  /* For each hemisphere S_{mu,lambda}, there is an associated element g of PSL_2(Z[w]), */
  /* which maps it onto a hemisphere of Floege's extended fundamental domain \hat{G}. */
  /* If g maps S_{mu,lambda} onto itself, we divide out this symmetry here. */

  local( Fraction, conj_mu_by_mu, RealPart );

  for( j= 1, numberOfSpheres,

    /* Exclude the case of the unit hemispheres, because they are in the corners */
    /* of the diagram, thus the cells are automatically divided. */
    if ( norm( nfbasistoalg(K,Mu[j])) > 1,

        /* Check the condition that Fraction := (lambda^2 +1)/mu is in the ring of */
        /* integers R. By Floege's lemma 10.4, this is equivalent to the existence of */
        /* a matrix g which mirrors the 2-cell with respect to a geodesic half-plane */
        /* passing through the hemisphere center. In this case, */
        /* g = [ -lambda, -(lambda^2 +1)/mu; mu, lambda].  */
        Fraction = nfeltdiv(K, nfeltmul(K, Lambda[j], Lambda[j])+[1,0]~, Mu[j]);         
        if( frac( Fraction) == [0,0]~,

          /* The matrix g reflects S_{mu,lambda} with respect to a half-plane */
          /* through the hemisphere center lambda/mu, */
          /* of constant real part, if conj(mu)/mu = 1 ; */
          /* of constant imaginary part, if conj_mu_by_mu := conj(mu)/mu = -1. */
          conj_mu_by_mu = nfeltdiv(K, conj( nfbasistoalg(K, Mu[j])), Mu[j]);
          if ( conj_mu_by_mu == [1,0]~,

            /* mirroring plane of constant real part */

            if( mIs3mod4,

              /* Check that the hemisphere center doesn't lie on the rim of the base rectangle. */
              RealPart = component(sphereCenter[j],1) -1/2*component(sphereCenter[j],2);
              if( RealPart > -1/2 && RealPart < 1/2,
                /* Exclude that the hemisphere center lies on an integer translate of */
                /* the imaginary axis, because the cell gets divided automatically. */
                if( frac( RealPart) != 0,

                        /* Create middle edge to cut off half of the 2-cell. */
                        insertVerticalEdge(j, RealPart);
                );
              );
            ,/* else m not congruent 3 modulo 4 */

              /* Check that the hemisphere center doesn't lie on the rim of the base rectangle. */
              RealPart = component(sphereCenter[j],1);
              if( RealPart > 0 && RealPart != 1/2 && RealPart < 1 ,
                /* Exclude that the hemisphere center lies on a 1/2 integer translate of */
                /* the imaginary axis, because the cell gets divided automatically. */
                if( frac( RealPart) != 0,

                        /* Create middle edge to cut off half of the 2-cell. */
                        insertVerticalEdge(j, RealPart);
                );
              );
            );
          ,/* else, mirroring plane of constant imaginary part */
            print("****Error in procedure prepareDivideOutAssociatedMatrices: horizontal mirroring needs to be implemented.");
          );
          if ( norm( nfbasistoalg(K, conj_mu_by_mu)) != 1, 
            error("conj_mu_by_mu not plusminus 1 as determined by Floege. Check procedure prepareDivideOutAssociatedMatrices");
          ); 
        );
    );
  );
  /* Renew the set of all vertices, as a union of the pointsOfSphere sets: */
  totalPointSet = Set(); 
  for (j = 1, numberOfSpheres,
        totalPointSet = setunion(totalPointSet, Set(pointsOfSphere[j]:list));
  );
};
 /* end of procedure prepareDivideOutAssociatedMatrices */


divideOutAssociatedMatrices() = 
{
  /* For each hemisphere S_{mu,lambda}, there is an associated element g of PSL_2(Z[w]), */
  /* which maps it onto a hemisphere of Floege's extended fundamental domain \hat{G}. */
  /* If g maps S_{mu,lambda} onto itself, we divide out this symmetry here. */

  local( Fraction, conj_mu_by_mu, RealPart );

  for( j= 1, numberOfSpheres,

    /* Exclude the case of the unit hemispheres, because they are in the corners */
    /* of the diagram, thus the cells are automatically divided. */
    if ( norm( nfbasistoalg(K,Mu[j])) > 1,

        /* Check the condition that Fraction := (lambda^2 +1)/mu is in the ring of */
        /* integers R. By Floege's lemma 10.4, this is equivalent to the existence of */
        /* a matrix g which mirrors the 2-cell with respect to a geodesic half-plane */
        /* passing through the hemisphere center. In this case, */
        /* g = [ -lambda, -(lambda^2 +1)/mu; mu, lambda].  */
        Fraction = nfeltdiv(K, nfeltmul(K, Lambda[j], Lambda[j])+[1,0]~, Mu[j]);         
        if( frac( Fraction) == [0,0]~,

          /* The matrix g reflects S_{mu,lambda} with respect to a half-plane */
          /* through the hemisphere center lambda/mu, */
          /* of constant real part, if conj(mu)/mu = 1 ; */
          /* of constant imaginary part, if conj_mu_by_mu := conj(mu)/mu = -1. */
          conj_mu_by_mu = nfeltdiv(K, conj( nfbasistoalg(K, Mu[j])), Mu[j]);
          if ( conj_mu_by_mu == [1,0]~,

            /* mirroring plane of constant real part */

            if( mIs3mod4,

              /* Check that the hemisphere center doesn't lie on the rim of the base rectangle. */
              RealPart = component(sphereCenter[j],1) -1/2*component(sphereCenter[j],2);
              if( RealPart > -1/2 && RealPart < 1/2,
                /* Exclude that the hemisphere center lies on an integer translate of */
                /* the imaginary axis, because the cell gets divided automatically. */
                if( frac( RealPart) != 0,

                        print("Vertically mirrored hemisphere at ", sphereCenter[j],
                                " of radius square ",   radiusSquare[j]);
                        addVertexIdentifications(j, Fraction);
                        verticalCellHalfDivision(j);
                );
              );
            ,/* else m not congruent 3 modulo 4 */

              /* Check that the hemisphere center doesn't lie on the rim of the base rectangle. */
              RealPart = component(sphereCenter[j],1);
              if( RealPart > 0 && RealPart != 1/2 && RealPart < 1,
                /* Exclude that the hemisphere center lies on a 1/2 integer translate of */
                /* the imaginary axis, because the cell gets divided automatically. */
                if( frac( RealPart) != 0,

                        print("Vertically mirrored hemisphere at ", sphereCenter[j],
                                " of radius square ",radiusSquare[j]);
                        addVertexIdentifications(j, Fraction);
                        verticalCellHalfDivision(j);
                );
              );
            );
          );
          if( conj_mu_by_mu == [-1,0]~,

                /* mirroring plane of constant imaginary part */ 
                print("Horizontally mirrored hemisphere at ", sphereCenter[j]," of radius square ",
                        radiusSquare[j]);
                addVertexIdentifications(j, Fraction);
                print("****Error in procedure divideOutAssociatedMatrices: horizontal cell subdivision needed.");
          );
          if ( norm( nfbasistoalg(K, conj_mu_by_mu)) != 1, 
                error("conj_mu_by_mu not plusminus 1 as determined by Floege. Check procedure divideOutAssociatedMatrices");
          ); 
        );
    );
  );
};
/* end of procedure divideOutAssociatedMatrices */



verticalCellHalfDivision(l) = 
{ /* In 2-cell number l, we cut off the corners of inferior real part as the sphere center. */

  local( RealPart, z, rsquare, cornerRealPart, j);
  local( positiveHalfCellCorners:list);
  positiveHalfCellCorners = listcreate( maxNeighboursNumber +2);

  j = twoCellSupport:list[l];

  /* Record the real part of the hemisphere center. */
  if( mIs3mod4,
      RealPart = component(sphereCenter[j],1) -1/2*component(sphereCenter[j],2);
  ,/* else */
      RealPart = component(sphereCenter[j],1);
  );
        
  for( k = 1, length( pointsOfSphere[j]:list),

        /* Consider the vertex (z,r) := pointsOfSphere[j][k] */
        z = component( pointsOfSphere[j]:list[k], 1);
        rsquare = component( pointsOfSphere[j]:list[k], 2);

        if( mIs3mod4,
          cornerRealPart = component(z,1) -1/2*component(z,2);
        , /* else */
          cornerRealPart = component(z,1);
        );
           
        if( cornerRealPart >= RealPart,
                /* Put it into the corners list of the positive half cell (to be kept). */ 
                listput( positiveHalfCellCorners, setsearch( totalPointSet,
                                                pointsOfSphere[j]:list[k]));
                print("point kept as corner at ",z," = ",nfbasistoalg(K,z) );
        );
  );
  cornersOf2cell[l] = positiveHalfCellCorners;
  listkill( positiveHalfCellCorners);
};
/* end of procedure verticalCellHalfDivision */



addPoint( pointAtIntersection, j, r) =
{
  local( RealPart, ImaginaryPartBySqrtm);

  RealPart = component(pointAtIntersection,1);
  ImaginaryPartBySqrtm = component( pointAtIntersection, 2);

  /* Check that pointAtIntersection is in the base rectangle: */
  if( 0 <= RealPart, 
        if( RealPart <= 1, 
          if( 0 <= ImaginaryPartBySqrtm, 
            if( ImaginaryPartBySqrtm <= 1, 
        
                /* Put pointAtIntersection into pointsOfSphere, pointsOfLine: */
                addVertex( pointAtIntersection, j);
                addVertexToLine( pointAtIntersection, j, length(linesOfSphere[j]));
                addVertexToLine( pointAtIntersection, j, r);
            );
          );
        );
  );
};
/* end of procedure addPoint */


addPoint3mod4( pointAtIntersection, j, r) =
{
  local( RealPart, ImaginaryPartBySqrtm);

  RealPart = component(pointAtIntersection,1) -1/2*component(pointAtIntersection,2);
  ImaginaryPartBySqrtm = 1/2*component( pointAtIntersection, 2);

  /* Check that pointAtIntersection is in the base rectangle: */
  if( -1/2 <= RealPart, 
        if( RealPart <= 1/2, 
          if( 0 <= ImaginaryPartBySqrtm, 
            if( ImaginaryPartBySqrtm <= 1/2, 
        
                /* Put pointAtIntersection into pointsOfSphere, pointsOfLine: */
                addVertex( pointAtIntersection, j);
                addVertexToLine( pointAtIntersection, j, length(linesOfSphere[j]));
                addVertexToLine( pointAtIntersection, j, r);
            );
          );
        );
  );
};
/* end of procedure addPoint3mod4 */




subdivideVerticalLoop() = 
{
  /* For m = 21, the vertical spurious edge of the two-cell lying on the sphere of radius 1/2, */
  /* is a loop: Its origin is equivalent to its end. */
  /* To be able to guarantee that also this cell is fixed pointwise by its stabilizer, */
  /* we subdivide it into two edges, if the above situation occurs. */
  /* We obtain this cellular subdivision by inserting a spurious horizontal edge. */ 
  local( RealPart);

  /* Pick the sphere of radius 1/2. */
  for( j = 1, numberOfSpheres,
     if( radiusSquare[j] == 1/4,

        /* Pick the line on which lies the the vertical spurious edge. */
        for( r = 1, length( linesOfSphere[j]),  
           /* Check if the actual line is vertical: */
           /* It suffices to check if the real part of the direction vector vanishes. */
           if( mIs3mod4,
                RealPart = component( component( linesOfSphere[j][r], 2), 1)
                      -1/2*component( component( linesOfSphere[j][r], 2), 2);
           , /* else not congruent to 3 mod 4 */
                RealPart = component( component( linesOfSphere[j][r], 2), 1);
           );
           if( RealPart == 0,

                if( IsLoop( pointsOfLine[j][r]:list),

                        performEdgeFission( j, r);
                );
           );
        );
     );
  );
  /* Renew the set of all vertices, as a union of the pointsOfSphere sets: */
  totalPointSet = Set(); 
  for (j = 1, numberOfSpheres,
        totalPointSet = setunion(totalPointSet, Set(pointsOfSphere[j]:list));
  );
};
/* end of procedure subdivideVerticalLoop */




performEdgeFission(j, r) = 
{ 
  /* We enter a Line:= linesOfSphere[j][r] with two points on it. */
  /* It is the underlying geometric structure of an edge with two vertices. */ 
  /* A copy of the Line is made. */
  /* One of the two points on the Line is moved to onto the copy. */
  /* The barycenter of the two points is put on both the Line and also on its copy. */
  /* We obtain the underlying geometric structure for two edges with a common vertex. */

  local( Line, auxiliaryVector, auxiliaryList, pointOne, pointTwo, Barycenter, heightsquare);

  Line = linesOfSphere[j][r];
  pointOne = pointsOfLine[j][r][1];
  pointTwo = pointsOfLine[j][r][2];
  Barycenter = ( component( pointOne, 1) +component( pointTwo,1) )/2;
  /* Compute the height of the hemisphere number j over Barycenter. */
  heightsquare = (1 -norm( nfbasistoalg(K, nfeltmul(K, Mu[j], Barycenter) -Lambda[j]))) 
                   /norm(nfbasistoalg(K,Mu[j]));

  /* Put a copy of the Line into linesOfSphere[j]. */
  auxiliaryList = linesOfSphere[j];
  listput(auxiliaryList,Line);
  linesOfSphere[j] = auxiliaryList;
  listkill( auxiliaryList);

  /* Insert a pointsOfLine list containing pointTwo and the Barycenter. */
  auxiliaryList = List([ pointTwo, [Barycenter, heightsquare] ]);
  if( length(pointsOfLine[j]) < length(linesOfSphere[j]),
                /* There are now more lines attributed to sphere number j, then */
                /* pointsOfLine lists. So we have to concatenate the new list. */
                auxiliaryVector = vector(1);
                auxiliaryVector[1] = auxiliaryList;
                auxiliaryVector = concat( pointsOfLine[j], auxiliaryVector);
                pointsOfLine[j] = auxiliaryVector;
  ,/* else there is enough space. */
                pointsOfLine[j][length(linesOfSphere[j])] = auxiliaryList;
  );

  /* The Line is now charged with pointOne and the Barycenter. */
  pointsOfLine[j][r] = List([ pointOne, [Barycenter, heightsquare] ]);

  /* Put the Barycenter into pointsOfSphere[j]. */
  auxiliaryList = pointsOfSphere[j]:list;
  listput( auxiliaryList, [Barycenter, heightsquare]);
  pointsOfSphere[j]:list = auxiliaryList;
};
/* end of procedure performEdgeFission */



insertVerticalEdge(j, RealPart) = 
{ /* A vertical edge is inserted on the line with constant real part RealPart */
  /* and intersected with the lines of sphere number j. */

  local( auxiliaryVector, auxiliaryList);

  /* Put the arc above the line with constant real part RealPart, */
  /* represented by a line through [RealPart, 0]~ of direction sqrt(-m) = x, */
  /* into linesOfSphere[j]. */
  auxiliaryList = linesOfSphere[j];
  auxiliaryList = concat(auxiliaryList, List([ [[RealPart,0]~, nfalgtobasis(K, x)] ]));

  linesOfSphere[j] = auxiliaryList;
  listkill( auxiliaryList);

  /* Insert an empty pointsOfLine list. */
  auxiliaryList = listcreate(maxNeighboursNumber+2);
  if( length(pointsOfLine[j]) < length(linesOfSphere[j]),
                /* There are now more lines attributed to sphere number j, then */
                /* pointsOfLine lists. So we have to concatenate the new empty list. */
                auxiliaryVector = vector(1);
                auxiliaryVector[1] = auxiliaryList;
                auxiliaryVector = concat( pointsOfLine[j], auxiliaryVector);
                pointsOfLine[j] = auxiliaryVector;
  ,/* else */
                pointsOfLine[j][length(linesOfSphere[j])] = auxiliaryList;
  );
  joinSpuriousVerticesOnHorizontalEdges(j, RealPart);
};
/* end of procedure insertVerticalEdge */




joinSpuriousVerticesOnHorizontalEdges(j, RealPart) =
{
  /* If not already in pointsOfSphere[j], we record a top, respectively a bottom vertex, */
  /* which lie on the intersection of */
  /* the boundary of the cell number j with its vertical fission axis. */
  local( pointAtIntersection, setOfPointsOfSphere);

  setOfPointsOfSphere = getSetOfPointsOfSphere(j);

  /* Intersect the new vertical line with all the horizontal lines in linesOfSphere[j]: */
  for (r = 1, length(linesOfSphere[j]) -1,

     /* Find out if this one is horizontal. */
     if( IsHorizontal( linesOfSphere[j][r]),

        pointAtIntersection = intersectLines( j, r, length(linesOfSphere[j]) );
        /* The vertical line is the last entry in linesOfSphere[j]. */

        /* Check that this vertex is not already in pointsOfSphere[j]. */
        if( setsearch( setOfPointsOfSphere, pointAtIntersection, 1),

                /* If so, record it into pointsOfSphere[j] and pointsOfLine[j][r]. */
                if( mIs3mod4,
        
                        addPoint3mod4( pointAtIntersection, j, r);
                ,/* else */             
                        addPoint( pointAtIntersection, j, r);
                );
        );
     );
  );
  recordPointsOntoAxis(j, RealPart);
};
/* end of procedure joinSpuriousVerticesOnHorizontalEdges */


recordPointsOntoAxis(j, RealPart) =
{
  local( examinedVertex, RealPartOfVertex);

  for( r = 1, length( pointsOfSphere[j]:list),

    examinedVertex = component( pointsOfSphere[j]:list[r], 1);
    if( mIs3mod4, 

        RealPartOfVertex = component(examinedVertex,1) -1/2*component(examinedVertex,2);

    , /* else m not congruent 3 mod 4 */

        RealPartOfVertex = component(examinedVertex,1);
    );
    if( RealPartOfVertex == RealPart,  

        /* Record it onto the vertical axis. */
        addVertexToLine( examinedVertex, j, length(linesOfSphere[j]));
        /* The vertical line is the last entry in linesOfSphere[j]. */
    );
  );
};
/* end of procedure recordPointsOntoAxis */



prepareSpuriousCellFissionNot3mod4() = 
{ /* This is a procedure for the NOT 3 modulo 4 case. */ 
  /* For each hemisphere with real part of its center equalling one half, */
  /* we prepare a fission along the vertical line of real part 1/2 into two cells. */

  local( RealPart);

  for( j= 1, numberOfSpheres,

      /* Check that the hemisphere center lies on the imaginary axis. */
      RealPart = component(sphereCenter[j],1);

      if( RealPart == 1/2,

        /* Put the arc above the imaginary axis, */
        /* represented by a vertical line through the origin c , */
        /* into linesOfSphere[j]. */
        insertVerticalEdge(j, RealPart);
      );
  );
  /* Renew the set of all vertices, as a union of the pointsOfSphere sets: */
  totalPointSet = Set(); 
  for (j = 1, numberOfSpheres,
        totalPointSet = setunion(totalPointSet, Set(pointsOfSphere[j]:list));
  );
};
/* end of procedure prepareSpuriousCellFissionNot3mod4 */


spuriousCellFissionNot3mod4() = 
{ /* This is a procedure for the NOT 3 modulo 4 case. */ 
  /* For each hemisphere of center with real part 1/2, */
  /* we perform a fission along the line of real part 1/2 into two cells. */
  /* We define the negative half cell as a new cell, */
  /* and keep only the non-negative corners and edges in the existing cell. */
  /* Already, an edge has been inserted along the line of real part 1/2 */
  /* by the procedure prepareSpuriousCellFissionNot3mod4(). */ 

  local( RealPart, z, rsquare);
  local( negativeHalfCellCorners:list, positiveHalfCellCorners:list);
  negativeHalfCellCorners = listcreate( maxNeighboursNumber +2);
  positiveHalfCellCorners = listcreate( maxNeighboursNumber +2);

  for( j= 1, numberOfSpheres,

      /* Check that the hemisphere center lies on the line of real part 1/2. */
      RealPart = component(sphereCenter[j],1);

      if( RealPart == 1/2,

        /* Record the corners with negative real part in negativeHalfCellCorners. */

        for( k = 1, length( pointsOfSphere[j]:list),

           /* Consider the vertex (z,r) := pointsOfSphere[j][k] */
           z = component( pointsOfSphere[j]:list[k], 1);
           rsquare = component( pointsOfSphere[j]:list[k], 2);

           if( component(z,1) <= 1/2,
                /* Put it into the corners list of the negative half cell. */ 
                listput( negativeHalfCellCorners, setsearch( totalPointSet,
                                                pointsOfSphere[j]:list[k]));
                /* print("negative point at ",z," = ",nfbasistoalg(K,z) ); */
           );
           if( component(z,1) >= 1/2,
                /* Put it into the corners list of the positive half cell. */ 
                listput( positiveHalfCellCorners, setsearch( totalPointSet,
                                                pointsOfSphere[j]:list[k]));
                /* print("positive point at ",z," = ",nfbasistoalg(K,z) ); */
           );
        );      

        createNegativeHalfCell( j, negativeHalfCellCorners:list);
        cornersOf2cell[j] = positiveHalfCellCorners;
        listkill( positiveHalfCellCorners);
        listkill( negativeHalfCellCorners);
     );
  );
};
/* end of procedure spuriousCellFissionNot3mod4() */



prepareSpuriousCellFission() = 
{ /* This is a procedure for the 3 modulo 4 case. */ 
  /* For each hemisphere centered on the imaginary axis, */
  /* we prepare a fission along this axis into two cells. */

  local( RealPart);

  for( j= 1, numberOfSpheres,

      /* Check that the hemisphere center lies on the imaginary axis. */
      RealPart = component(sphereCenter[j],1) -1/2*component(sphereCenter[j],2);

      if( RealPart == 0,

        /* Put the arc above the imaginary axis, */
        /* represented by a vertical line through the origin c , into linesOfSphere[j]. */
        insertVerticalEdge(j, RealPart);
      );
  );
  /* Renew the set of all vertices, as a union of the pointsOfSphere sets: */
  totalPointSet = Set(); 
  for (j = 1, numberOfSpheres,
        totalPointSet = setunion(totalPointSet, Set(pointsOfSphere[j]:list));
  );
};
/* end of procedure prepareSpuriousCellFission */


spuriousCellFission() = 
{ /* This is a procedure for the 3 modulo 4 case. */ 
  /* For each hemisphere centered on the imaginary axis, */
  /* we perform a fission along this axis into two cells. */
  /* We define the negative half cell as a new cell, */
  /* and keep only the non-negative corners and edges in the existing cell. */
  /* Already, an edge has been inserted along the imaginary axis */
  /* by the procedure prepareSpuriousCellFission(). */ 

  local( RealPart, z);
  local( negativeHalfCellCorners:list, positiveHalfCellCorners:list);
  negativeHalfCellCorners = listcreate( maxNeighboursNumber +2);
  positiveHalfCellCorners = listcreate( maxNeighboursNumber +2);

  for( j= 1, numberOfSpheres,

      /* Check that the hemisphere center lies on the imaginary axis. */
      RealPart = component(sphereCenter[j],1) -1/2*component(sphereCenter[j],2);

      if( RealPart == 0,

        /* Record the corners with negative real part in negativeHalfCellCorners. */

        for( k = 1, length( pointsOfSphere[j]:list),

           /* Consider the vertex (z,r) := pointsOfSphere[j][k] */
           z = component( pointsOfSphere[j]:list[k], 1);

           if( component(z,1) -1/2*component(z,2) <= 0,
                /* Put it into the corners list of the negative half cell. */ 
                listput( negativeHalfCellCorners, setsearch( totalPointSet,pointsOfSphere[j]:list[k]));
                /* print("negative point at ",z," = ",nfbasistoalg(K,z) ); */
           );
           if( component(z,1) -1/2*component(z,2) >= 0,
                /* Put it into the corners list of the positive half cell. */ 
                listput( positiveHalfCellCorners, setsearch( totalPointSet,
                                                pointsOfSphere[j]:list[k]));
                /* print("positive point at ",z," = ",nfbasistoalg(K,z) ); */
           );
        );      

        createNegativeHalfCell( j, negativeHalfCellCorners:list);
        cornersOf2cell[j] = positiveHalfCellCorners;
        listkill( positiveHalfCellCorners);
        listkill( negativeHalfCellCorners);
     );
  );
};
/* end of procedure spuriousCellFission */




createNegativeHalfCell( j, negativeHalfCellCorners:list) =
{
  numberOf2cells++;
  listput( twoCellSupport:list, j);
  cornersOf2cell[ numberOf2cells] = negativeHalfCellCorners;
};
/* end of procedure createNegativeHalfCell */




getVertexOrbits() =
{
  local( counter:small, entrysum:small, areIdentified:small);

  vertexOrbitNumber = vector( length(totalPointSet));
  vertexOrbitRepresentative = vector( ceil( length(totalPointSet)/2)+3);
  counter = 1;

  for( k = 1, length(totalPointSet),
    entrysum = 0;

    for( l = 1, k-1, /* run through the points which have already been assigned an orbit number. */
        areIdentified = listSearch(equivalentVertices[k],l);

        if( entrysum == 0 && areIdentified > 0,
                /* Put the point number k into the orbit of point number l. */
                vertexOrbitNumber[k] = vertexOrbitNumber[l];

                chooseVertexRepresentative(k,l);
        );
        entrysum = entrysum + areIdentified;
    );
    if ( entrysum == 0, /* point number k is the first to occur in its orbit */
        vertexOrbitRepresentative[counter] = k;
        vertexOrbitNumber[k] = counter;
        counter = counter +1;
    );
  );
  logfile = Str( logfile,"There are ",counter-1," orbits of vertices.");
  numberOfVertexOrbits = counter -1;
};
/* end of the procedure getVertexOrbits. */



chooseVertexRepresentative(k,l) =
{
  /* Choose for the orbit representative the vertex amongst */
  /* the k'th one and the existing representative, */
  /* the one whose projection onto C of is nearer to zero; */
  /* provided that it is in the preferred stripe of width 1/2. */
  local( oldRepresentative, point_k, inStripe);
                
  oldRepresentative = component(eval(totalPointSet[
                                vertexOrbitRepresentative[vertexOrbitNumber[l]]]),1);
  point_k = component(eval(totalPointSet[k]),1);

  inStripe = checkWhichPointsAreInStripe( oldRepresentative, point_k);
                
  if( inStripe == "k",
        vertexOrbitRepresentative[vertexOrbitNumber[l]] = k;
  );
  if( inStripe == "both" || inStripe == "none",

        /* Choose for the orbit representative the vertex */
        /* whose projection is nearer to the origin. */

        if( norm(nfbasistoalg(K, oldRepresentative)) > norm(nfbasistoalg(K, point_k)),

                /* Take the actual point as orbit representative. */
                vertexOrbitRepresentative[vertexOrbitNumber[l]] = k;
        );
  );
  /* else, in the case "old", keep the old representative. */
};
/* end of procedure chooseVertexRepresentative */



chooseEdgeRepresentative(k,l) =
{
  /* Choose for the orbit representative the edge amongst */
  /* the k'th one and the existing representative, */
  /* the barycenter of whose projection onto C of is nearer to zero; */
  /* provided that the latter one is in the preferred stripe of width 1/2. */
  local( oldBarycenter, Barycenter_k, inStripe);
                
  oldBarycenter = (component(eval(component(
                        edgesList:list[edgeOrbitRepresentative[ edgeOrbitNumber[l]]],1)),1)
                 +component(eval(component(
                        edgesList:list[edgeOrbitRepresentative[ edgeOrbitNumber[l]]],2)),1) 
                )/2;
  Barycenter_k = ( component(eval(component(edgesList:list[k],1)),1)
                  +component(eval(component(edgesList:list[k],2)),1) )/2;

  inStripe = checkWhichPointsAreInStripe( oldBarycenter, Barycenter_k);

  if( inStripe == "k",
        edgeOrbitRepresentative[edgeOrbitNumber[l]] = k;
  );
  if( inStripe == "both" || inStripe == "none",

        if( norm(nfbasistoalg(K, oldBarycenter)) > norm(nfbasistoalg(K, Barycenter_k)),

                        edgeOrbitRepresentative[edgeOrbitNumber[l]] = k;
        );
  );
  /* else, in the case "old", keep the old representative. */
};
/* end of procedure chooseEdgeRepresentative */





getEdgeOrbits() =
{
  local( counter:small, entrysum:small, areIdentified:small);
  /* Use the vector of lists equivalentEdges to sort the edges into their orbits. */

  edgeOrbitNumber = vector( length(edgesList:list)+2+m);
  edgeOrbitRepresentative = vector( length(edgesList:list)+1+m);
  counter = 1;

  for( k = 1, length(edgesList:list),
    entrysum = 0;

    for( l = 1, k-1,

        areIdentified =  listSearch(equivalentEdges[k],l) 
                        +listSearch(equivalentEdges[l],k);

        if( entrysum == 0  && areIdentified > 0,
                /* Put the edge number k into the orbit of edge number l */
                edgeOrbitNumber[k] = edgeOrbitNumber[l];

                /* Choose for the orbit representative the edge amongst */
                /* the k'th one and the existing representative, */
                /* the barycenter of whose projection onto C of is nearer to zero; */
                /* provided that the latter one is in the preferred stripe of width 1/2. */
                chooseEdgeRepresentative(k,l);
        );
        entrysum = entrysum +areIdentified ;
    );
    if ( entrysum == 0, /* edge number k is the first to occur in its orbit */
        edgeOrbitRepresentative[counter] = k;
        edgeOrbitNumber[k] = counter;
        counter++;
    );
  );
  logfile = Str( logfile,"There are ",counter -1," orbits of edges.");
  numberOfEdgeOrbits = counter -1;
};
/* end of procedure getEdgeOrbits */



sortByOrbitNumbers() = 
{
  /* Sort the entries in the vectors cornersOf2cell by orbit numbers. */
  local( orbitNumbersVector, numberOfCorners, newCornersOfSphere);

  for ( j = 1, numberOf2cells,
        numberOfCorners = length(cornersOf2cell[j]);
        orbitNumbersVector = vector(numberOfCorners);
        newCornersOfSphere = vector(numberOfCorners);

        /* Record the orbit numbers of the corners in cornersOf2cell[j]: */
        for ( r = 1, numberOfCorners,
                orbitNumbersVector[r] = vertexOrbitNumber[cornersOf2cell[j][r]];
        );

        orbitNumbersVector = vecsort(orbitNumbersVector,,1);
         /* Binary digit 1 of flag means: 
          indirect sorting, return the permutation instead of the permuted vector. */

        /* Constitute a copy of cornersOf2cell[j], but with the ordering given
                by the permutation vector: */
        for ( r = 1, numberOfCorners,
                newCornersOfSphere[r] = cornersOf2cell[j][orbitNumbersVector[r]];
        );
        cornersOf2cell[j] = newCornersOfSphere;
  );
};
/* end of procedure sortByOrbitNumbers */


swapCorners(k,s,t) =
{
  /* Swap the corners s and t in cornersOf2cell[k]. */
  local( auxiliaryVector);
  auxiliaryVector = cornersOf2cell[k];
  auxiliaryVector[t] = cornersOf2cell[k][s];
  auxiliaryVector[s] = cornersOf2cell[k][t];
  cornersOf2cell[k] = auxiliaryVector;
};
/* end of procedure swapCorners */



checkTwoCellIdentification(j,k) = 
{
  local( barycentricallyIdentified, jointlyIdentified);

  /* Check if the two 2-cells number j and k are identified by isometries from Gamma. */
  jointlyIdentified = checkIdentificationIntersection(j,k);
  barycentricallyIdentified = checkBarycenterIdentification(j,k);

  if(  jointlyIdentified == 0 && barycentricallyIdentified,
        print("***Error detected in 2-cell structure: 2-cells number ",j," and ",k,
                " barycentrically but not jointly identified.");
  );
  if(  barycentricallyIdentified == 0 && jointlyIdentified,
        print("***Error detected in 2-cell structure: 2-cells number ",j," and ",k,
                " jointly but not barycentrically identified.");
  );

  if ( jointlyIdentified && barycentricallyIdentified,

        /* If so, delete one of them. */
        chooseRepresentative2cellAmongst(j,k);  
  );
};
 /* end of procedure checkTwoCellIdentification */



chooseRepresentative2cellAmongst(j,k) =
{
  /* Delete the 2cell amongst numbers k and j (on the same orbit), whose */
  /* barycenter is on the negative half-plane, has real part greater than 1/2, */
  /* or has the bigger distance to the origin. */
  local( inStripe, barycenter_j, barycenter_k );

  barycenter_j = getBarycenter(j);
  barycenter_k = getBarycenter(k);

  inStripe = checkWhichPointsAreInStripe( barycenter_j, barycenter_k);

  if( inStripe == "k",

        deleteCellFlag[j] = 1;  
  );    
  if( inStripe == "old",

        deleteCellFlag[k] = 1;          
  );
  if( inStripe == "both" || inStripe == "none",

     if( mIs3mod4,

                chooseCellNearerToAxis(j,k, barycenter_j,  barycenter_k);

     ,/* else  not congruent 3 mod 4 */  
        /* Delete the 2-cell amongst numbers k and j (on the same orbit), whose */
        /* barycenter has the bigger distance to the origin. */
        if ( norm( nfbasistoalg(K, barycenter_j )) > norm( nfbasistoalg(K, barycenter_k )),
  
                deleteCellFlag[j] = 1; 
                , /* else */
                deleteCellFlag[k] = 1;
        );
     );
  );
}; /* end of procedure chooseRepresentative2cellAmongst(j,k) */



chooseCellNearerToAxis(j,k, barycenter_j,  barycenter_k) =
{
  local( RealPart_k, RealPart_j);

  RealPart_k  = component( barycenter_k, 1)  -component( barycenter_k, 2)/2;
  RealPart_j  = component( barycenter_j, 1)  -component( barycenter_j, 2)/2;

  if( abs(RealPart_k) >   abs(RealPart_j),
 
          deleteCellFlag[k] = 1; 
  , /* else */
        if( abs(RealPart_k) <   abs(RealPart_j),

                deleteCellFlag[j] = 1;
        , /* else  they have equal real part, and we delete the one farther to the origin.*/
                if ( norm( nfbasistoalg(K, barycenter_j )) 
                   > norm( nfbasistoalg(K, barycenter_k )),
  
                        deleteCellFlag[j] = 1; 
                        , /* else */
                        deleteCellFlag[k] = 1;
                );
        );
  );
}; /* end of procedure chooseCellNearerToAxis */



initialize2cells()=
{
  /* To start, put a 2-cell onto each kept hemisphere. */
  /* Later, we will pass to the quotient by the action, and keep only one 2-cell per orbit. */
  twoCellSupport = listcreate( numberOfSpheres +m);
  for( j = 1, numberOfSpheres,
        listput( twoCellSupport:list, j);
  );
  deleteCellFlag = vector(numberOfSpheres +m);
  numberOf2cells = numberOfSpheres;
  cornersOf2cell = vector(numberOfSpheres +m);
};
/* end of procedure initialize2cells */



recordCorners() =
{ 
  /* Get the totalPointSet-numbers of the points in pointsOfSphere 
                and record them into cornersOf2cell. */
  for (j = 1, numberOfSpheres,
    cornersOf2cell[j] = vector(length(pointsOfSphere[j]:list));
    for( r = 1, length(pointsOfSphere[j]:list),
        cornersOf2cell[j][r] = setsearch( totalPointSet, 
                                pointsOfSphere[j]:list[r]);
    );
  );
};
/* end of procedure recordCorners */


identify2cells(orbitSum) =
{  
  /* Check which two-cells are identified by the action of Gamma. */
  for (j = 1, numberOf2cells,
    for (k = j+1, numberOf2cells,
      /* Check if the vertices on the hemispheres number j and k are on the same orbits:*/
      if( orbitSum[j] == orbitSum[k],
        checkTwoCellIdentification(j,k);
      );
    );
  );
};
/* end of procedure identify2cells */



putIntoedgesOf2cell(j, newEdgeNumber) =
{
        /* put the new edge into edgesOf2cell[j] */
        local(auxiliaryList);
        auxiliaryList = edgesOf2cell[j];
        listsort( auxiliaryList, 1); /* cleanse list of double occurencies. */
        listput(auxiliaryList, newEdgeNumber);
        edgesOf2cell[j] = auxiliaryList; 
};
/* end of procedure putIntoedgesOf2cell */


getOriginAndEndNumbers() =
{ 
  /* Record the numbers of the origin and end points of the edges */
  local( representativeOriginNumber);

  EdgeOrigin = vector( length(edgesList:list));
  EdgeEnd = vector( length(edgesList:list));

  /* For every edge in edgesList, record the following information: */  
  for ( j = 1,  length(edgesList:list),

        /* record number in totalPointSet of the origin vertex */
        EdgeOrigin[j] = setsearch( totalPointSet, component( edgesList:list[j], 1));

        /* Check if this is origin is on the orbit of the origin of the edge's representative. */
        representativeOriginNumber = setsearch( totalPointSet, component( edgesList:list[ 
                                        edgeOrbitRepresentative[ edgeOrbitNumber[j]]], 1));
        if( listSearch(equivalentVertices[ EdgeOrigin[j]], representativeOriginNumber) > 0,

                /* record number in totalPointSet of the end vertex */
                EdgeEnd[j] = setsearch( totalPointSet, component( edgesList:list[j], 2));

        , /* else we must invert the edge to give the induced orientation to it. */

                EdgeEnd[j] = EdgeOrigin[j];
                EdgeOrigin[j] = setsearch( totalPointSet, component( edgesList:list[j], 2));
        );
  );
}; 
/* end of procedure getOriginAndEndNumbers */