\\ ---------------  GP code  ------------------------------------------------------------
\\ 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.6 of June 30, 2011.    
\\---------------------------------------------------------------------------------------

global( logfile, K, x, m, mIs3mod4, w, numberOfLines:small, maxNeighboursNumber:small, totalPointSet, pointsOfSphere, pointsOfLine, edgesList:list=List([]), twoCellSupport:list=List([]), sphereCenter, radiusSquare, neighbours, linesOfSphere, stabilizer, edgeStabilizer, identificationMatrices, equivalentVertices, transportFrom, equivalentEdges, edgesOf2cell, EdgeOrigin, EdgeEnd, vertexOrbitNumber, vertexOrbitRepresentative, numberOfVertexOrbits:small, numberOfEdgeOrbits:small, edgeOrbitNumber, edgeOrbitRepresentative, Mu:list=List([]),Lambda:list=List([]), numberOfSpheres:small, cornersOf2cell, deleteFlag:list=List([]), numberOf2cells, deleteCellFlag, boundaryMatrix, boundaryMap, numberOfTwoCells, CellBoundaryMatrix, singularpoint:list=List([]), singularDenominator:list=List([]), classNumber, EdgeTwoDim, edge2Label:list, EdgeThreeDim, edge3Label:list); 



print("\n Loading the script Bianchi.gp :");
print("Program for the computation of the quotient of hyperbolic three-space");
print("by the Bianchi groups PSL_2(R), R a ring of imaginary quadratic integers.")
print("\n Copyright (C) 2010 by Alexander D. Rahm. \n Bianchi.gp is a free software covered by the GNU General Public License.");
print("Bianchi.gp version 2.0.6 of June 30, 2011. \n");
print("Please use Pari/GP version 2.4.3 or more recent:");
if (lex(version(), [2,4,3]) < 0, error("Please use Pari/GP version 2.4.3 or more recent."););
print("\n The program is started with the command Start(m), \n","where m has to be chosen amongst the positive square-free integers.");
print("The imaginary quadratic number field K is then defined as the field of the");
print("rational numbers extended by x with x^2 = -m, and R is its ring of integers.");



\r Bianchi-common_procedures_and_diagnosis.gp
\r Bianchi-Fundamental_Polyhedron_functions.gp
\r Bianchi-Fundamental_Polyhedron_procedures.gp
\r Bianchi-orbit_space_functions.gp
\r Bianchi-orbit_space_procedures.gp
\r Bianchi-group_homology.gp
\r Bianchi-output.gp


Start(mInput) =   /* User call of the main program. */{
  m = mInput;
  if( issquarefree(m), 

    MainProgram();
    print("\n","To save the computed data in backup files in the current folder, type the command BackUp()");
    print("To output a selection of the obtained information in .tex files, type the command OutputFiles()"); 
  ,/* else */
        print(m,"=m is not squarefree. ",
            "We would get a subgroup of the Bianchi group for m/(squares in m).");
  );            
}; /* end of the program Bianchi.gp. */





fdomain(mInput) = /* Call the main program and backup. */{
  m = mInput;
  if (issquarefree(m), 

    MainProgram();
    if ( m >= 10,  /* Create backup files for the computed data. */
                BackUp( "record");
        );
 ,/* else */
        print(m,"=m is not squarefree. ",
            "We would get a subgroup of the Bianchi group for m/(squares in m).")
  );            
}; /* end of the program Bianchi.gp. */




MainProgram() =
{
  gettime(); 
  logfile = " ";
  initializeNumberField(m);
  getSingularPoints();
  getFundamentalDomain();

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

        compareHeights(3);
        createLines();
        getVertices(); 
        prepareSpuriousCellFissionNot3mod4();
        compareHeights(0);
        removePoints(); 
        removeLines(); 
        logfile = Str( logfile,"m=",m," Cell structure found in ",gettime()," milliseconds.");
        if( abs(K.disc) > 400, BackUp("fundamental domain"););
        ComputeOrbitSpaceNot3mod4();

  ,/* else m congruent 3 mod 4 */
        compareHeights(1);
        createLines();
        getVertices();
        removePoints(); 
        removeLines(); 
        compareHeights(2);
        createLines();
        getVertices();
        compareHeights(3);
        createLines();  
        getVertices();
        prepareSpuriousCellFission();
        compareHeights(0);
        prepareDivideOutAssociatedMatrices();
        removePoints();
        removeLines();
        logfile = Str( logfile, "m=",m," Cell structure found in ",gettime()," milliseconds.");
        if( abs(K.disc) > 400, BackUp("fundamental domain"););
        ComputeOrbitSpace3mod4();
  );  
  logfile = Str( logfile,"m=",m," orbit space found in ",gettime()," milliseconds.");
  logfile = Str( logfile, "m=",m," gives ", length(edgesList)," edges and ", 
                        length(totalPointSet), " vertices in the base rectangle.");
  if( m!= 1 && m != 3,
    ComputeGroupHomology();
  ); 
 }; 
/* end of the main program of Bianchi.gp. */
 





initializeNumberField( mInput=m) =
{
  m = mInput; 

  if (m<0, error("m negative. Want to use Q[sqrt(-m)]"));

  /* Define the number field K globally. */
  K = nfinit( x^2 +m);
  w = nfbasistoalg(K, [0,1]~);
  mIs3mod4 = (Mod(m,4) == Mod(3,4));
  classNumber = qfbclassno(K.disc);
  
  print("Starting the computation of the Floege PSL_2(Z[w])-cellular complex, with w = ",w,".");      

  if ( mIs3mod4,
        if( w != Mod(1/2*x - 1/2, x^2 + m),
                error(Str("Base should be chosen such that w = ", Mod(1/2*x -1/2, x^2+ m),
                        ". Please use version 2.4.3 or higher of Pari/GP for the case m =",m,"."));
        );
  );
}; /* end of the procedure initializeNumberField. */



BackUp( mode="record") =
{
  local( namesOfVariables);
  namesOfVariables = [ "numberOfLines", 
        " maxNeighboursNumber", "totalPointSet", "pointsOfSphere", 
        " pointsOfLine", " edgesList", " sphereCenter", 
        "radiusSquare", " neighbours", " linesOfSphere", 
         " stabilizer", "edgeStabilizer", 
        "identificationMatrices", "equivalentVertices", "transportFrom", "equivalentEdges", 
        "edgesOf2cell", " EdgeOrigin","EdgeEnd", 
        "vertexOrbitNumber", " vertexOrbitRepresentative", 
        "numberOfVertexOrbits", "numberOfEdgeOrbits",
        " edgeOrbitNumber", "edgeOrbitRepresentative", "Mu", "Lambda", 
        " numberOfSpheres", "numberOf2cells", "twoCellSupport",
        "cornersOf2cell", "deleteFlag", "deleteCellFlag",
        "boundaryMatrix", " boundaryMap", " numberOfTwoCells",
        "CellBoundaryMatrix", "singularpoint", " singularDenominator"];

  if( mode == "record",
        system(Str("mkdir Bianchi-",m));
        for( j = 1, length( namesOfVariables),
                write( Str( "Bianchi-",m,"/m=", m, namesOfVariables[j], ".txt"),
                        eval(namesOfVariables[j]) 
                );
        );
  );
  if( mode == "fundamental domain",
        system(Str("mkdir fundamental_domain-",m));
        for( j = 1, length( namesOfVariables),
                write( Str( "fundamental_domain-",m,"/m=", m, namesOfVariables[j], ".txt"),
                        eval(namesOfVariables[j]) 
                );
        );
  );
  if( type(mode) == "t_INT",

    m = mode;
    initializeNumberField(m);
    getSingularPoints();

    for( j = 1, length( namesOfVariables),

        /*Attention, do not try to read new variables which are not yet recorded. */      
        eval( Str( namesOfVariables[j], "=", 
                   read( Str( "m=", m, namesOfVariables[j], ".txt")) 
                  )
             );
    );
    print("You can now use the command ComputeGroupHomology()");
  );
};
addhelp(BackUp, "To load data for the case m which is stored in the current folder, type BackUp(m)"); 
/* end of procedure BackUp. */


Restore(mInput) = {
  ReadData(mInput);
  convertIdentificationMatrices();
  computeRemainingOrbitSpaceStructure();
};
addhelp(Restore, "To load old data (Bianchi.gp-1.x format) for the case m which is stored in the current folder, type Restore(m)"); 




ReadData(mInput) = {
  local( namesOfVariables);
  namesOfVariables = [ "numberOfLines", 
        " maxNeighboursNumber", "totalPointSet", "pointsOfSphere", 
        " pointsOfLine", " edgesList", " sphereCenter", 
        "radiusSquare", " neighbours", " linesOfSphere", 
         " stabilizer", "edgeStabilizer", 
        "vertexIdentifications", " edgesOf2cell", " EdgeOrigin","EdgeEnd",
        "vertexOrbitNumber", " vertexOrbitRepresentative", 
        "numberOfVertexOrbits", "numberOfEdgeOrbits",
        " edgeOrbitNumber", "edgeOrbitRepresentative", "Mu", "Lambda", 
        " numberOfSpheres", "numberOf2cells", "twoCellSupport",
        "cornersOf2cell", "deleteFlag", "deleteCellFlag",
        "boundaryMatrix", " boundaryMap", " numberOfTwoCells",
        "CellBoundaryMatrix", "singularpoint", " singularDenominator"];

  m = mInput;
  initializeNumberField(m);
  getSingularPoints();

  for( j = 1, length( namesOfVariables),

        /*Attention, do not try to read new variables which are not yet recorded. */      
        eval( Str( namesOfVariables[j], "=", 
                   read( Str( "m=", m, namesOfVariables[j], ".txt")) 
                  )
             );
  );
};


convertIdentificationMatrices() = {
  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, vertexIdentifications[r,r] );
      transportFrom[r] = sublistPut(transportFrom[r], length(identificationMatrices));
      equivalentVertices[r] = sublistPut( equivalentVertices[r], r);    
  );
  for( r = 1, length(totalPointSet),

    for( s = r+1, length(totalPointSet),
        
        if( length(vertexIdentifications[r,s]) > 0,

            listput( identificationMatrices, vertexIdentifications[r,s] );
            transportFrom[r] = sublistPut(transportFrom[r], length(identificationMatrices));
            equivalentVertices[r] = sublistPut( equivalentVertices[r], s);

            listput( identificationMatrices, vertexIdentifications[s,r] );
            transportFrom[s] = sublistPut(transportFrom[s], length(identificationMatrices));
            equivalentVertices[s] = sublistPut( equivalentVertices[s], r);
        );
    );
  );
    /* Free memory from old format entries: */ vertexIdentifications = 0;
};


computeRemainingOrbitSpaceStructure() = {
  recordSingularSelfIdentifications(recordVectorOfIdentifications()); 
  for(j = 1, 4,
    recordSingularTransitivityRelations();
  );
  if( mIs3mod4,
        getVertexOrbits();
        initialize2cells();
        recordCorners();
        divideOutAssociatedMatrices();
        spuriousCellFission();
        sortByOrbitNumbers(); 
        orbitSum = recordOrbitSum();
        getEdges(); 
        computeEdgeStabilizers();
        getEdgeIdentifications();
        getEdgeOrbits();
        identify2cells( orbitSum);
        getOriginAndEndNumbers();
        computeEdgeStabilizers();

  , /* else m is not congruent to 3 mod 4 */
        initialize2cells();
        getVertexOrbits();
        recordCorners();
        spuriousCellFissionNot3mod4();
        divideOutAssociatedMatrices();
        sortByOrbitNumbers();
        orbitSum = recordOrbitSum();
        identify2cells( orbitSum);
        getEdges();
        computeEdgeStabilizers();
        getEdgeIdentifications();
        getEdgeOrbits();
        getOriginAndEndNumbers();
        computeEdgeStabilizers();
  );
  print("You can now use the command ComputeGroupHomology()");
};