\\ ---------------  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.1 of June 28th, 2010.    
\\---------------------------------------------------------------------------------------



getSingularPoints() =
{
  local( sSquare);
  print("The following singular points are in the boundary of ",
                "the Bianchi fundamental polyhedron: \n");
/* The singular points are given by p*(r+x)/s, where p runs through
   the respective function "threeModFour" or "oneOrTwoModFour" below. */
/* The latter functions are implementations of a formula by R. G. Swan. */

  singularpoint       = listcreate(m*(classNumber-1));
  singularDenominator = listcreate(m*(classNumber-1));

  for (s = 2, m,
        sSquare = s^2;
        for (r = floor(1 -(s/2)), floor(s/2),
                if (sSquare <= r^2 +m, 
                        if ( mIs3mod4,
                                threeModFour(r,s)
                                ,
                                oneOrTwoModFour(r,s)
                        );
                );      
        );
  );
  if ( length(singularpoint) == 0, 
        print("There are no singular points.");
       ,
        listsort(singularpoint);
        listsort(singularDenominator);
        if ( length(singularpoint) != length(singularDenominator),
                error("Ordering singular points list failed in procedure getSingularPoints.");
        );
  );
  print();
};  
{addhelp(getSingularPoints, "Compute singular points of imaginary quadratic field K = Q[x] with -x^2 = m not congruent to 0 mod 4. Call this function with argument a positive integer m.");}

threeModFour(r,s) =
{
 local( point, equivalentPoint, sHalf);
 if ( s > 2,
  if ( Mod(s,2) == Mod(0,2),
    sHalf = s/2;
    if ( Mod(r^2 +m, 2*s) == Mod(0, 2*s),
        for (p = 1, sHalf -1,
                if ( gcd(p, sHalf) == 1,

                        point = orbitrepresentative(p,r,s);
                        listput(singularpoint, point);
                        print("p=",p," r=",r," s=",s, 
                                " give a singular point at ", point);
                        listput(singularDenominator, [point, s] );

                        if (component( point, 1) == 0,

                          equivalentPoint = nfalgtobasis(K, 1 +p*(r+x)/s);
                          print("which has another orbit representative at ",equivalentPoint);
                          listput(singularpoint, equivalentPoint);
                          listput(singularDenominator, [equivalentPoint, s] );
                        );

                        if (component(point, 1) >= 1/2,

                          equivalentPoint = point -[1,0]~;
                          print("which has another orbit representative at ",equivalentPoint);
                          listput(singularpoint, equivalentPoint);
                          listput(singularDenominator, [equivalentPoint, s] );
                        );

                        if (component(point, 1) -component( point, 2)/2 == 1/2,

                          equivalentPoint = nfalgtobasis(K, p*(r+x)/s -1);
                          print("which has another orbit representative at ",equivalentPoint);
                          listput(singularpoint, equivalentPoint);
                          listput(singularDenominator, [equivalentPoint, s] );
                        );
                );
        );
    );
  );
 );
};

oneOrTwoModFour(r,s) =
{
   local( point, equivalentPoint);

   if ( Mod(r^2 +m, s) == Mod(0, s),
        for (p = 1, s-1,
                if ( gcd(p,s) == 1,

                        point = orbitrepresentative(p,r,s);
                        listput(singularpoint, point);
                        print("p=",p," r=",r," s=",s," give a singular point at ",point);
                        listput(singularDenominator, [point, s] );

                        if ( component(point, 1) == 0,

                          equivalentPoint = nfalgtobasis(K, 1+p*(r+x)/s);
                          print("which has another orbit representative at ",equivalentPoint);
                          listput(singularpoint, equivalentPoint);
                          listput(singularDenominator, [equivalentPoint, s] );
                        );
                );
        );
   );
};

orbitrepresentative(p,r,s) =
{
  local(PR);

  PR = lift( Mod( p*r, s));

  /*Return */  nfalgtobasis( K, (PR +p*x)/s)
}; 
/* end of the procedure getSingularPoints(). */



getBarycenter(j) =
{
  local(barycenter);
  /* Construct the barycenter of the projections of the vertices on the 2-cell number j onto the complex plane. */
  barycenter = 0;
  for ( r = 1, length( cornersOf2cell[j]),
        barycenter = barycenter 
                        + component( eval( totalPointSet[cornersOf2cell[j][r]]), 1);
  );
  barycenter = barycenter /  length( cornersOf2cell[j]);

  /* return */ barycenter
};
/* end of function getBarycenter */


createPreliminaryMu(mufixabssquare) =
{               
  local( aSquare, b,  squareRoot, preliminaryMu:list);  

  preliminaryMu = listcreate(4000);

  /* Create the list preliminaryMu of denominators Mu of cusps; Mu will be [a,b]~,
  such that the norm square of Mu = a^2 +m*b^2 = mufixabssquare.*/
        
  for ( a = 0, floor(sqrt(mufixabssquare)),
        /* We do not run through negative values for a, because */
        /* for a = Re(Mu) < 0, it suffices to consider -Mu. */
    
          aSquare = a^2;
          squareRoot = sqrt( (mufixabssquare -aSquare)/m );

          b = round( +squareRoot); /* positive root */

          if ( aSquare +m*b^2 == mufixabssquare,  
            
                listput(preliminaryMu, [a,b]~);
          );
          b = round( -squareRoot); /* negative root */

          if ( aSquare +m*b^2 == mufixabssquare,  
            
                listput(preliminaryMu, [a,b]~);
          );
  );
  /*Return the list  preliminaryMu. */ preliminaryMu
};
  {addhelp(createPreliminaryMu,"Create list preliminaryMu of denominators Mu of hemisphere centers; Mu will be [a,b]~. Call this function with argument the square of the norm of Mu.");}
/* end of the function createPreliminaryMu(mufixabssquare). */




createPreliminaryMu3mod4(mufixabssquare) =
{               
  local( aSquare, b, squareRoot, preliminaryMu:list);   

  preliminaryMu = listcreate(4000);

  /* Create the list preliminaryMu of denominators Mu of cusps; Mu will be [a,b]~,
  such that the norm square of Mu = a^2 -a*b +b^2 *(1+m)/4 = mufixabssquare.*/

  for ( a = 0, floor(sqrt(mufixabssquare)),
        /* We do not run through negative values for a, because */
        /* for a = Re(Mu) < 0, it suffices to consider -Mu. */
    
        aSquare = a^2;
        squareRoot = sqrt( (1+m)*mufixabssquare -m*aSquare);

        b = round( 2/(1+m) *(a +squareRoot)); /* +sqrt(Discriminant) */
                                         
        if ( mufixabssquare == aSquare -a*b +b^2 *(1+m)/4, /* == normsquare([a,b]~) */ 

                listput(preliminaryMu, [a,b]~);
        );
        b = round( 2/(1+m) *(a -squareRoot)); /* -sqrt(Discriminant) */
                                         
        if ( mufixabssquare == aSquare -a*b +b^2 *(1+m)/4,  /* == normsquare([a,b]~) */

                listput(preliminaryMu, [a,b]~);
        );
  );
  /*Return the list  preliminaryMu. */ preliminaryMu
};
/* end of the function createPreliminaryMu3mod4(mufixabssquare). */




addVertex( pointAtIntersection, j) =
{ /* Write pointAtIntersection into the lists pointsOfSphere[j]. */

  local(jzetasquare, auxiliaryList:list);

  jzetasquare = (1 -norm( nfbasistoalg(K, nfeltmul(K, Mu[j], pointAtIntersection) -Lambda[j])))
                /norm(nfbasistoalg(K,Mu[j]));
  /* jZeta is the height of the hemisphere j */
  /* (and of the one at number neighbours[j][k] for the suitable k) over pointAtIntersection. */

  if (jzetasquare >= 0,
                        
        /* Write pointAtIntersection into pointsOfSphere[j]: */
        auxiliaryList = pointsOfSphere[j]:list;
        
        /* Check that pointAtIntersection is not already in pointsOfSphere[j]: */ 
        if( setsearch(Set(auxiliaryList) ,[pointAtIntersection,jzetasquare], 1),

                listput(auxiliaryList, [pointAtIntersection, jzetasquare]);
        );
        listsort(auxiliaryList,1); 

        /* Delete points which are below another point in auxiliaryList. */
        if( length(auxiliaryList) > 1,
            for  (r = 1, length(auxiliaryList)-1,

                if (component(auxiliaryList[r],1) == component(auxiliaryList[r+1],1),

                        if(component(auxiliaryList[r],2) > component(auxiliaryList[r+1],2),

                                     error("Ordering auxiliaryList failed in procedure addVertex");
                        );
                        listput(auxiliaryList, auxiliaryList[r+1],r);
                );
            );
        );
        pointsOfSphere[j] = auxiliaryList;
        listkill(auxiliaryList);
  ); 
};
/* end of procedure addVertex */



addVertexToLine( pointAtIntersection, j, k) =
{ /*Write pointAtIntersection into the lists pointsOfLine[j][k]. */

  local(jzetasquare, auxiliaryList:list);

  jzetasquare = (1 -norm( nfbasistoalg(K, nfeltmul(K, Mu[j], pointAtIntersection) -Lambda[j])))
                /norm(nfbasistoalg(K,Mu[j]));
  /* jZeta is the height of the hemispheres j and neighbours[j][k] over pointAtIntersection. */

  if (jzetasquare >= 0,
                        
        /*Write pointAtIntersection into pointsOfLine[j][k]: */
        auxiliaryList = pointsOfLine[j][k]:list;
        listput(auxiliaryList, [pointAtIntersection, jzetasquare]);
        listsort(auxiliaryList,1);

        /* Delete points which are below another point in auxiliaryList. */      
        if( length(auxiliaryList) > 1,
            for  (s = 1, length(auxiliaryList)-1,

                if (component(auxiliaryList[s],1) == component(auxiliaryList[s+1],1),

                        if(component(auxiliaryList[s],2) > component(auxiliaryList[s+1],2),

                                     error("auxiliaryList disorder");
                        );
                        listput(auxiliaryList, auxiliaryList[s+1],s);
                );
            );
        );
        pointsOfLine[j][k] = auxiliaryList;
        listkill(auxiliaryList);
  ); 
};


algMatrix( givenMatrix) = { 
  /* Convert the entries of the givenMatrix, which are given in a basis */
  /* of the ring of integers, into the "POLMOD" format. This is the */
  /* algebraic notation, as a polynomial in x modulo x^2 +m. */
  /* This allows us to use Pari/GP matrix multiplication with */
  /* the converted matrix. */
  local( convertedMatrix); convertedMatrix = matrix(2,2);

  for ( j = 1, 2,
        for( k = 1, 2,
                convertedMatrix[j,k] = nfbasistoalg(K, givenMatrix[j,k]);
        );
  );
  /* return */ convertedMatrix
}; /* end of the function algMatrix. */


algMatrixList( matrixList:list) = {
  /* Convert a list of matrices, with entries given in the basis */
  /* of the ring of integers, into a set of matrices with */ 
  /* entries in the "POLMOD" format. This is the */
  /* algebraic notation, as a polynomial in x modulo x^2 +m. */
  local( convertedList: list);

  convertedList = listcreate( length( matrixList));

  for(j = 1, length( matrixList),

        listput( convertedList, algMatrix( matrixList[j]) );
  );
  /* return */ convertedList
}; /* end of the function algMatrixList. */


matrixToBasis( givenMatrix) = { 
  /* Convert the entries of the givenMatrix, which are given in a basis */
  /* of the ring of integers, into the "POLMOD" format. This is the */
  /* algebraic notation, as a polynomial in x modulo x^2 +m. */
  /* This allows us to use Pari/GP matrix multiplication with */
  /* the converted matrix. */
  local( convertedMatrix); convertedMatrix = matrix(2,2);

  for ( j = 1, 2,
        for( k = 1, 2,
                convertedMatrix[j,k] = nfalgtobasis(K, givenMatrix[j,k]);
        );
  );
  /* return */ convertedMatrix
}; /* end of the function algMatrix. */


matrixListToBasis( matrixList:list) = {
  /* Convert a list of matrices, with entries given in the basis */
  /* of the ring of integers, into a set of matrices with */ 
  /* entries in the "POLMOD" format. This is the */
  /* algebraic notation, as a polynomial in x modulo x^2 +m. */
  local( convertedList: list);

  convertedList = listcreate( length( matrixList));

  for(j = 1, length( matrixList),

        listput( convertedList, matrixToBasis( matrixList[j]) );
  );
  /* return */ convertedList
}; /* end of the function algMatrixList. */


setToList( incomingSet) = {
  local(outputList:list);
  outputList = listcreate(length(incomingSet));

  for( j = 1, length(incomingSet),

        listput(outputList, eval(incomingSet[j]));
  );
  /* return */ outputList
};


Composition(firstMatrixSet, secondMatrixSet) = {
  local( firstMatrices:list, secondMatrices:list, composedMatrices:list);
  firstMatrices = algMatrixList(setToList(firstMatrixSet));
  secondMatrices = algMatrixList(setToList(secondMatrixSet));
  composedMatrices = listcreate(min(10,length(firstMatrices)*length(secondMatrices)));

  for( j = 1, min(5,length(firstMatrices)),
        for( k = 1, min(5,length(secondMatrices)),
        
           if( length(composedMatrices) < 10,
                listput( composedMatrices, firstMatrices[j]*secondMatrices[k]);
           );
        );
  );
  /*  suppress all but one occurence of each element in the list */
  composedMatrices = matrixListToBasis( setToList(Set(composedMatrices)));
  /* return */ Set(composedMatrices)
};


\\ ---------------  GP code  ------------------------------------------------------------
\\ Diagnosis 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 1.0.0 of October 29th, 2010.   
\\---------------------------------------------------------------------------------------




checkFor2cellInversions() = /* inactive check, to be called from command line */
{
  /* The matrices which rotate a 2-cell onto itself should all have been found */
  /* by Floege's algorithm. This procedure, which is not part of the program, */
  /* gives an indication if this has indeed been achieved. */
  /* We display the vertex orbits of the corners of a 2-cell as soon as there */
  /* is one which occurs more than once. */

  local( numberOfOrbits, orbitsList:list);

  for(  j = 1, numberOf2cells,

        /* record the number of distinct vertex orbits on which are the corners */
        /* of 2-cell number j. */
        orbitsList = listcreate( length( cornersOf2cell[j]));

        for( r = 1, length( cornersOf2cell[j]),

                listput( orbitsList, vertexOrbitNumber[ cornersOf2cell[j][r]]);
        );
        listsort( orbitsList, 1);
        numberOfOrbits = length( orbitsList);
        /* the flag "1" in listsort is for deleting double entries */

        if( numberOfOrbits < length( cornersOf2cell[j]),

                print("The following vertex orbits occur for corners of sphere number ",j,": ",
                        orbitsList );
        );
  );
};
/* end of diagnosis procedure checkFor2cellInversions */



checkFor2cellstabilizers() = /* inactive check, to be called from command line */
{
  /* The stabilizers of the 2-cells should be trivial. To check this, */
  /* we search for 2cells whose edge stabilizers all have the same cardinality. */

  local( edgeStabilizerCardinals:list, numberOfCardinals);

  for(  j = 1, numberOf2cells,

        /* record the number of distinct vertex orbits on which are the corners */
        /* of 2-cell number j. */
        edgeStabilizerCardinals = listcreate( length( edgesOf2cell[j]));
        for( r = 1, length( edgesOf2cell[j]),

                listput( edgeStabilizerCardinals, length( edgeStabilizer[ edgesOf2cell[j][r]]));
        );
        numberOfCardinals = length( listsort( edgeStabilizerCardinals, 1));
        /* the flag "1" in listsort is for deleting double entries */

        if( numberOfCardinals == 1 && edgeStabilizerCardinals[1] > 2,
                /* allow for the trivial stabilizer +/-1 of cardinal 2 */
                print("The following cardinal occurs for all edge stabilizers of sphere number "
                        ,j,": ", edgeStabilizerCardinals[1] );
        );
  );
};
/* end of diagnosis procedure checkFor2cellstabilizers */


checkPointwiseFixing() = 
{  
  /* Perform a check if all the cells in the cellular complex are fixed pointwise */
  /* by their stabilizer. This is obviously true for vertices. */
  /* We want to check if some edge can be sent onto itself reverting its orientation. */
  /* This is possible when origin and end of the edge lie on the same orbit. */
  /* In this case, we produce an error message calling for a subdivision of the edge. */

  for ( r = 1, length( edgesList),
        if( vertexOrbitNumber[ EdgeOrigin[r]] == vertexOrbitNumber[ EdgeEnd[r]],
                print("*** Error detected by procedure checkPointwiseFixing: Edge number ",r," should be subdivided.");
        );
  );
  checkPointwiseFixing2cells();
};
/* end of diagnosis procedure checkPointwiseFixing. */

checkPointwiseFixing2cells() = 
{
  /* The matrices which rotate a 2-cell onto itself should all have been found */
  /* by Floege's algorithm. This procedure checks that this has indeed been achieved. */
  /* We make sure that each 2-cell has at least three corners whose vertex orbits */
  /* do not appear twice for this 2-cell. */
  /* Only two such uniquely occuring representatives could lie on a symmetry axis. */

  local( numberOfUniquelyOccuringOrbits, numberOfCorners, orbitsList:list);

  for(  j = 1, numberOf2cells,

        numberOfUniquelyOccuringOrbits = 0;

        /* Record the vertex orbits on which lie the corners of 2-cell number j. */
        orbitsList = listcreate( length( cornersOf2cell[j]));
        for( r = 1, length( cornersOf2cell[j]),

                listput( orbitsList, vertexOrbitNumber[ cornersOf2cell[j][r]]);
        );
        listsort( orbitsList);
        numberOfCorners = length( orbitsList);
        if( numberOfCorners != length( cornersOf2cell[j]) || numberOfCorners < 3,
                print("***Error detected by procedure checkPointwiseFixing2cells: incorrect number of corners on 2-cell ",j);
        );
        if( orbitsList[1] != orbitsList[2],
                numberOfUniquelyOccuringOrbits++;
        );
        for( r = 2, numberOfCorners-1,
                if( orbitsList[r] != orbitsList[r-1] && orbitsList[r] != orbitsList[r+1],
                        numberOfUniquelyOccuringOrbits++;
                );
        );
        if( orbitsList[numberOfCorners-1] != orbitsList[numberOfCorners],
                numberOfUniquelyOccuringOrbits++;
        );
        if( numberOfUniquelyOccuringOrbits < 3,
                print("***Warning:  2-cell ",j," not guaranteed by procedure checkPointwiseFixing2cells to be pointwise fixed.");
        );
        listkill( orbitsList);
  );
};
/* end of diagnosis procedure checkPointwiseFixing2cells. */



checkIdealClass_cusp_correspondence() =
{
  /* Check the correspondence between the ideal classes and the cusps in the quotient of */
  /* Hyperbolic 3-space by the action of the Bianchi group, by counting singular cusps. */
  local( numberOfSingularPoints, heightSquare);

  /* Count the singular cusps in the quotient space. */
  numberOfSingularPoints = 0;
  for ( j = 1, numberOfVertexOrbits,
        heightSquare = component( eval(totalPointSet[ vertexOrbitRepresentative[j]]), 2);
        if( heightSquare == 0,
                numberOfSingularPoints ++;
        );
  );
  if( numberOfSingularPoints +1 == classNumber,
        print("There are ", numberOfSingularPoints," singular cusps in the quotient space, ",
                "representing the non-trivial ideal classes.");
  ,/* else */
        print("***Error: There are ", numberOfSingularPoints," singular cusps in the quotient, ",
                "but ",classNumber -1, " non-trivial ideal classes.");
  );
};
/* end of diagnosis function checkIdealClass_cusp_correspondence */


findErrorInDifferential() =
{  
  local(A,B);
  for( j = 1, numberOf2cells, 
                
        A =  vecextract(CellBoundaryMatrix,Str(j)); 

        for( k = 1, numberOfVertexOrbits, 

                B = vecextract(boundaryMatrix ,Str(k),Str("1..",numberOfEdgeOrbits )); 

                if( B*A != 0, 

                        print("Differential ",B*A," at cell ",j,", vertex ",k);
                );
        ); 
  );
};
/* end of diagnosis function findErrorInDifferential */