\\ ---------------  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.9 of July 20, 2011.   
\\---------------------------------------------------------------------------------------



deleteSpheresBelow( preliminaryMu:list, preliminaryLambda:list, MuNumber:list, everywhere_below:list) =
{
  /* Remove preliminary hemispheres which are everywhere below
        some hemisphere already in the list. */
  local(MU:list,LAMBDA:list, sphere_center, invertedNormSquareOfMU, distanceSquare, a,b,c,j);
  MU = listcreate(numberOfSpheres);
  LAMBDA = listcreate(numberOfSpheres);
  sphere_center    = vector( length( preliminaryLambda));
  invertedNormSquareOfMU = vector( length( preliminaryLambda));

  print("Removing everywhere below ones amongst ",
        length(preliminaryLambda)," hemispheres...");

  /* Precompute sphere center and radius. */
  for (k:small = 1, length(preliminaryLambda),

    if (everywhere_below[k] == 0, 

        sphere_center[k] = nfeltdiv(K, preliminaryLambda[k],preliminaryMu[MuNumber[k]]);

        invertedNormSquareOfMU[k] = 1/norm(nfbasistoalg(K,preliminaryMu[MuNumber[k]]));                 
    );
  );

  /* Perform the inner loop of Algorithm 2. */
  for (k:small = 1, length(preliminaryLambda),

    j = 1;
    while ( everywhere_below[k] == 0 && j <= length(Lambda),

      if ( deleteFlag[j] == 0,

        distanceSquare = norm(nfbasistoalg(K,sphere_center[k] -sphereCenter[j]));

        /* Check if invertedNormOfMU[k] +distance <= radius[j].  Let */
        a = invertedNormSquareOfMU[k]; b = distanceSquare; c = radiusSquare[j];
        /* The condition sqrt(a) +sqrt(b) < sqrt(c) is equivalent to */
        /* 4*a*b < (c -a -b)^2 && a +b < c, because a, b, c > 0. */
        if( 4*a*b < (c -a -b)^2 && a +b < c, 

                        everywhere_below[k] = 1;
        );
        j++;
      );
    );
  );    
  listkill(MU); listkill(LAMBDA); listkill(preliminaryMu); 
  listkill(preliminaryLambda); listkill(MuNumber);
  /* return */ everywhere_below:list
}; 
{addhelp(deleteSpheresBelow,"Remove preliminary hemispheres which are everywhere below some hemisphere already in the list.")};
/* end of the function deleteSpheresBelow */


computeLine(j,k) =
{
  local(A,v, PositionVector,DirectionVector, Return, centerdifference);

  A = matrix(2,2);
  centerdifference = sphereCenter[neighbours[j][k]] -sphereCenter[j];

  if (mIs3mod4,
  
        A[1,1] = 2*component( centerdifference, 1) -component( centerdifference,2);

        A[1,2] = (m+1)/2 *component( centerdifference, 2) -component( centerdifference,1);

  
  , /* Else m not congruent 3 mod 4: */
        A[1,1] = 2*component( centerdifference,   1);
        A[1,2] = m*2*component( centerdifference, 2);
  );
  v = norm(nfbasistoalg(K,sphereCenter[neighbours[j][k]])) 
        - norm(nfbasistoalg(K,sphereCenter[j]))
        -radiusSquare[neighbours[j][k]] 
        + radiusSquare[j];

  if( A[1,1] == 0 && A[1,2] == 0,

        if ( v == 0,
                print("Accident: The two hemispheres ",j," and ",
                        neighbours[j][k], " are identic!");
                ,/* else */
                print("One hemisphere completely covers the other one");                
        );
  ,/* else */
        numberOfLines = numberOfLines +1;
        if (A[1,2] == 0,

                PositionVector = [v/A[1,1],0]~ ;
                DirectionVector = [0,1]~;

                /* line equation x= v/A[1,1] */

        ,/* Else, not a vertical line, as A[1,2] is nonzero: */

                /* line equation y = v/A[1,2] -(A[1,1]/A[1,2])*x */
                PositionVector = [0,v/A[1,2]]~ ;
                DirectionVector = [1,-A[1,1]/A[1,2]]~;  
        );
  );
  Return = [PositionVector, DirectionVector];
  Return
}; /* end of the function computeLine. */



intersectLines(j:small,k:small,r:small) =
{
  local( determinant, xi, direction_k, position_k, direction_r, position_r);

  direction_k = component( linesOfSphere[j]:list[k], 2);
  position_k  = component( linesOfSphere[j]:list[k], 1);
  direction_r = component( linesOfSphere[j]:list[r], 2);
  position_r  = component( linesOfSphere[j]:list[r], 1);
  xi = 0;

  determinant = -component(direction_k,1) *component(direction_r,2)
                +component(direction_k,2) *component(direction_r,1);

  if (determinant != 0,
      /* If the two lines have different slope, compute
                 pointAtIntersection.  */
        
        xi = 1/determinant * (-component(direction_r,2)
                *component(position_r-position_k,1)
                +component(direction_r,1)
                *component(position_r-position_k,2) );

  , /* else the two lines are parallel and have no intersection point. */
        /* We then return a vector which is surely outside the base rectangle. */
        position_k = [99,99]~;
  );
  /* Return pointAtIntersection = */ position_k+xi*direction_k
};
/* end of function intersectLines */


compareSides(printedgesonscreen=0) =
{
  local( edge, edgesCheckList1:list, edgesCheckList2:list, cancelCriterion);

  edgesCheckList1 = listcreate(2*numberOfLines);
  edgesCheckList2 = listcreate(2*numberOfLines);

  for (j = 1, numberOfSpheres,
   if( deleteFlag[j] == 0,

     if( printedgesonscreen == 1,       
             print("Cell ",j," lies on the hemisphere with center at ",
                  nfeltdiv(K,Lambda[j],Mu[j]) );
             print(" with radius square ", radiusSquare[j],
                " and has the following edges:");
      ); 

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

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

         if(length(edge) == 2, 

                if(printedgesonscreen == 1, print(edge); );

                if( setsearch(Set(edgesCheckList1),edge),

                        if(setsearch(Set(edgesCheckList2),edge),
                        
                                if (printedgesonscreen == 1,
                                  print("***Error in function compareSides: triple edge in cell diagram");
                                );
                               ,
                                listput(edgesCheckList2,edge);
                        );
                       , /* else not yet entered in edgesCheckList1 */
                        listput(edgesCheckList1,edge);
                );
                        
         );
         if(length(edge) > 2,
                print("***Error in function compareSides: edge with three corners");
         );
      );
   );   
  );

  if(Set(edgesCheckList1) == Set(edgesCheckList2),
     
        print("All ", length(edgesCheckList1), " = ",
                length(edgesCheckList2)," edges appear twice.");
        if( length(edgesCheckList1) > 3,
                cancelCriterion = 1
                ,
                cancelCriterion = 0;
        );
        , /* else some edges don't appear twice. */
                print("Some edges don't appear twice:");
                print(setminus(Set(edgesCheckList1),Set(edgesCheckList2)));
  );
  listkill(edgesCheckList1); listkill(edgesCheckList2);
  /* Return */ cancelCriterion
};
{addhelp(compareSides,"Prints the edge data if You call this as compareSides(1).");}
/* end of function compareSides. */



findUnitSphere(Center) =
{
  local( indexOfUnitSphere);
  /* Find the unit hemisphere with center "Center". */
  indexOfUnitSphere = 1;
  while ( sphereCenter[indexOfUnitSphere] != Center,
        indexOfUnitSphere++;
        if( indexOfUnitSphere > numberOfSpheres, 
                error("indexOfUnitSphere out of range. Compare vector sphereCenter and function findUnitSphere");
        );
  );
  /* Return the number of the found unit hemisphere: */   indexOfUnitSphere
};
/* end of function findUnitSphere */



estimateLimit() =
{
  /* Estimate the inferior limit for the hemisphere radius, */
  /* which is the reciprocal of the square of */
  /* the height of the lowest proper vertex of the fundamental domain. */
  /* If this estimation is not precise, it will be corrected in the below while loop. */
  /* Then the program will take more computation time. */
  local(limit);
  limit = vector(1879);
  limit[1]   = 1;
  limit[2]   = 4;
  limit[3]   = 1;
  limit[5]   = 25;
  limit[6]   = 34;
  limit[7]   = 1;
  limit[10]  = 66;
  limit[11]  = 6;    
  limit[13]  = 86;
  limit[14]  = 108;
  limit[15]  = 17;
  limit[17]  = 134;
  limit[19]  = 10;
  limit[21]  = 192;
  limit[22]  = 226;
  limit[23]  = 18;
  limit[26]  = 300;
  limit[29]  = 342;
  limit[30]  = 386;
  limit[31]  = 30;
  limit[32]  = 177;
  limit[33]  = 1091;
  limit[34]  = 482;
  limit[35]  = 41;
  limit[36]  = 226;
  limit[37]  = 534;
  limit[38]  = 588;
  limit[39]  = 49;
  limit[41]  = 646; 
  limit[43]  = 22;
  limit[47]  = 49;
  limit[51]  = 64;
  limit[55]  = 86;
  limit[58]  = 1282;
  limit[59]  = 39;
  limit[67]  = 34;
  limit[68]  = 134;
  limit[71]  = 93;
  limit[79]  = 95;
  limit[83]  = 49;
  limit[87]  = 121;
  limit[91]  = 100;
  limit[95]  = 121;
  limit[103] = 85;
  limit[107] = 121;
  limit[111] = 177;
  limit[115] = 144;
  limit[123] = 209;
  limit[139] = 103;
  limit[159] = 495;
  limit[163] = 82;
  limit[167] = 739;
  limit[179] = 161;
  limit[187] = 209;
  limit[193] = 9842; /* tentative */
  limit[211] = 121;
  limit[231] = 617;
  limit[235] = 324;
  limit[251] = 300;
  limit[267] = 801;
  limit[283] = 208;
  limit[307] = 207;
  limit[331] = 190;
  limit[379] = 217;
  limit[403] = 484;
  limit[427] = 576;
  limit[499] = 468;
  limit[547] = 488;
  limit[571] = 1200;
  limit[643] = 341;
  limit[667] = 801;
  limit[715] = 2081;
  limit[763] = 1409;
  limit[883] = 689;
  limit[907] = 704;
  limit[1027] = 1296;
  limit[1051] = 830;
  limit[1507] = 2409;
  limit[1879] = 3661;
  
  if( m <= 1879 && limit[m] > 0,

        limit = limit[m]; /* Pick the value based on computational experience, */

  , /* else use an estimation formula, based on extrapolation of the above values: */

        if ( mIs3mod4 == 0, /* m not congruent to 3 mod 4 */
                /* discriminant = -4m */ 
                limit = ceil(3*sqrt(2*classNumber)*4*m);
        , /* else m congruent 3 mod 4 */
                /* discriminant = -m */ 
                limit = ceil(3*sqrt(classNumber/2)*m);
        );
  );
  /* Return the estimated limit for m: */ limit
};
/* end of function estimateLimit */


Get_minimal_proper_vertex_heightsquare() =
{
  local( minzetasquare, z, zetasquare);

  minzetasquare = 1;

  for (j = 1, length(totalPointSet),
             
        z = component(eval(totalPointSet[j]),1);
        zetasquare = component(eval(totalPointSet[j]),2);
                
        if ( /* z not singular */ setsearch(Set(singularpoint),z,1),

                if ( minzetasquare > zetasquare,

                        minzetasquare = zetasquare;
                );
                if (zetasquare == 0,

                        print("***Error in function Get_minimal_proper_vertex_heightsquare: Hideous singular vertex at ", eval(totalPointSet[j]));
                );
        );
  );
  /* Return */ minzetasquare
};
/* end of function Get_minimal_proper_vertex_heightsquare(). */



increaseWithinNormsOfIntegers( mufixsquare) =
{
  local(F); /* Construct Set(F) to contain the norms of the elements of the ring of integers. */
  if( mIs3mod4,

        F = vecsort( concat( Vec( matrix( 2*round(sqrt(m))+3, m+3, 
                                        n, j, (j-1)^2 -(j-1)*(n-1) +((n-1)^2)*(1+m)/4 ))));
  , /* else */

        F = vecsort( concat( Vec( matrix( round(sqrt(m))+3, m+3, n, j, ((n-1)^2)*m+(j-1)^2))));
  );
  mufixsquare++;
        
  /* Increase mufixsquare until it meets an element of Set(F). */
  while( setsearch( Set(F), mufixsquare) == 0,

        mufixsquare++;
        if ( mufixsquare > vecmax(F),
                     error("mufixsquare out of range. Function increaseWithinNormsOfIntegers needs to be adapted.");
        );
  );
  /* Return */ mufixsquare
};
/* end of function increaseWithinNormsOfIntegers .*/