\\ \\ 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.1.1 of October 4, 2011. \\--------------------------------------------------------------------------------------- getFundamentalDomain() = { /* This procedure is the core of the program. */ local(mufixsquare, oldlength, limit, cancelCriterion, minzetasquare); system(Str("mkdir fundamental_domain-",m)); write( Str( "fundamental_domain-",m,"/Start_m", m,, ".txt"), "Starting the computation."); limit = estimateLimit(); /* Estimate the inferior limit for the hemisphere radius. */ /* If this estimation is not precise, it will be corrected in the below while loop. */ /* Then the program will take more computation time. */ Mu = listcreate(4*m*limit +16); Lambda = listcreate(4*m*limit +16); deleteFlag = listcreate (4*m*limit +16); oldlength = 0; cancelCriterion = 0; mufixsquare = 1; /* The value mufixsquare will control the radius of the Euclidean hemispheres. */ print("Running until the inverse of the radius square is greater than ", limit); while( cancelCriterion == 0, while( mufixsquare <= limit, print("Generating hemispheres of radius square = 1/",mufixsquare, " ..."); if ( mIs3mod4, createSpheres3mod4( mufixsquare); , createSpheres( mufixsquare); ); mufixsquare = increaseWithinNormsOfIntegers( mufixsquare); ); if (numberOfSpheres != oldlength, removeEverywhereBelowOnes(); cleanseListOfSpheres(); getLinesAndVertices(); minzetasquare = Get_minimal_proper_vertex_heightsquare(); oldlength = numberOfSpheres; , /* else numberOfSpheres == oldlength */ print("...kept old ",numberOfSpheres," hemispheres, ", length(totalPointSet)," vertices,"); ); if( minzetasquare*mufixsquare >= 1, print("All remaining vertices have height square at least ", minzetasquare,", "); print("and thus lie above the radius square of the ignored hemispheres, " "which is at most 1/",mufixsquare); cancelCriterion = 1; , /* else minzetasquare becomes the new expected lowest value for 1/mufixsquare. */ limit = 1/minzetasquare; print("Height denominator square difference between ignored hemispheres ", ,"and these vertices: "); print(mufixsquare," - ",1.0*limit); limit = min(limit, mufixsquare*3/2); ); ); system(Str("mkdir Bianchi-",m)); }; /* end of the procedure getFundamentalDomain(). */ cleanseListOfSpheres(keepPoints = 0) = /* Rewrite the list of hemispheres, omitting deleted entries. */ { local(MU:list, LAMBDA:list, j); MU = listcreate(numberOfSpheres); LAMBDA = listcreate(numberOfSpheres); for (k = 1, length(Lambda), if(deleteFlag[k] == 0, listput(MU, Mu[k]); listput(LAMBDA, Lambda[k]); ); ); if( keepPoints == 1, /* Shift the entries of the pointsOfSphere vector to their new hemisphere number. */ j = 1; for (k= 1, length(MU), while( deleteFlag[j] == 1, j++; ); pointsOfSphere[k]:list = pointsOfSphere[j]:list; j++; ); for( k= length(MU)+1, length(Lambda), pointsOfSphere[k]:list = List(); ); ); listkill(Mu); listkill(Lambda); listkill(deleteFlag); for (k= 1, length(MU), listput(Mu,MU[k]); listput(Lambda, LAMBDA[k]); listput(deleteFlag, 0); ); numberOfSpheres = length(Lambda); precomputeSphereData(); computeNeighbours(); listkill(MU); listkill(LAMBDA); }; {addhelp(cleanseListOfSpheres,"Rewrite the list of hemispheres, omitting deleted entries.");} /* end of the procedure cleanseListOfSpheres. */ precomputeSphereData() = /* Compute hemisphere center, radius, and list of neighbouring lines & hemispheres. */ { sphereCenter = vector(numberOfSpheres); radiusSquare = vector(numberOfSpheres); for( j = 1, numberOfSpheres, sphereCenter[j] = nfeltdiv(K,Lambda[j],Mu[j]); radiusSquare[j] = 1/norm(nfbasistoalg(K,Mu[j])); ); }; /* end of the procedure precomputeSphereData(). */ computeNeighbours() = { /* Generate the lists of neighbouring hemispheres */ local( a, b, c); neighbours = vector(numberOfSpheres); /* (the global object which is computed) */ for( j = 1, numberOfSpheres, neighbours[j] = vector(0); ); for( j = 1, numberOfSpheres, for( k = j+1, numberOfSpheres, /* Check that |sphereCenter[j] -sphereCenter[k]|^2 <= (radius[j] + radius[k])^2 */ c = norm(nfbasistoalg(K, sphereCenter[j] -sphereCenter[k])); a = radiusSquare[j]; b = radiusSquare[k]; /* The condition sqrt(c) <= sqrt(a) +sqrt(b) is equivalent to */ /* (a+b > c) || (a +b <= c && 4*a*b >= (c -a -b)^2), because a, b, c >= 0. */ if( (a+b > c) || (4*a*b >= (c -a -b)^2), neighbours[j] = concat(neighbours[j],[k]); neighbours[k] = concat(neighbours[k],[j]); ); ); ); /* And now compute the global maxNeighboursNumber. */ maxNeighboursNumber = 0; for( j = 1, numberOfSpheres, if ( length(neighbours[j]) > maxNeighboursNumber, maxNeighboursNumber = length(neighbours[j]); ); ); }; /* end of the procedure computeNeighbours(). */ createSpheres( mufixsquare) = { /* Version for the case m is not congruent 3 mod 4 */ local(preliminaryMu:list, preliminaryLambda:list, MU); local( j,k,numberOfMu:list,localDeleteFlag:list,nfFraction, Capacity, normOfMuPlus_wMu); Capacity = mufixsquare^2+100; preliminaryLambda = listcreate(Capacity); localDeleteFlag = listcreate (Capacity); numberOfMu = listcreate (Capacity); j=1; k=0; numberOfSpheres = 0; preliminaryMu = createPreliminaryMu( mufixsquare); while (j <= length(preliminaryMu) , MU = nfbasistoalg(K, preliminaryMu[j]); normOfMuPlus_wMu = round(sqrt(norm(MU +nfbasistoalg(K,nfeltmul(K,w,preliminaryMu[j]))))); /* Create Lambda := [a,b]~, which will be stored in the list preliminaryLambda. */ for ( a = -normOfMuPlus_wMu, normOfMuPlus_wMu, for ( b = -normOfMuPlus_wMu, normOfMuPlus_wMu, /* Check if Lambda/Mu is in the base rectangle: */ nfFraction = nfeltdiv(K, [a,b]~, preliminaryMu[j]); if( 0 <= component( nfFraction, 1), if( component( nfFraction, 1) <= 1, if( 0 <= component( nfFraction, 2), if( component( nfFraction, 2) <= 1, if ( idealadd(K, MU, [a,b]~) == [1,0;0,1], k = k+1; numberOfSpheres = numberOfSpheres +1; listput( preliminaryLambda, [a,b]~, numberOfSpheres); listput(numberOfMu,j,numberOfSpheres); listput(localDeleteFlag,0,numberOfSpheres); ); ); ); ); ); ); ); j++; k=0; ); if (length(Lambda) == 0, for (a= 1, length(preliminaryLambda), listput(Lambda,preliminaryLambda[a]); listput(Mu,preliminaryMu[numberOfMu[a]]); listput(deleteFlag,0); ); numberOfSpheres = length(Lambda); removeEverywhereBelowOnes(); cleanseListOfSpheres(); ); localDeleteFlag = deleteSpheresBelow( preliminaryMu,preliminaryLambda,numberOfMu,localDeleteFlag); /* write the global lists Lambda, Mu */ for (a = 1, length(preliminaryLambda), if (localDeleteFlag[a] == 0, listput(Lambda,preliminaryLambda[a]); listput(Mu,preliminaryMu[numberOfMu[a]]); listput(deleteFlag,0); ); ); numberOfSpheres = length(Lambda); listkill(preliminaryMu); listkill(preliminaryLambda); listkill(numberOfMu); listkill(localDeleteFlag); precomputeSphereData(); }; {addhelp(createSpheres,"Create the list Lambda with Lambda/Mu hemisphere center.");} /* end of the procedure createSpheres. */ createSpheres3mod4( mufixsquare) = { local(preliminaryMu:list, preliminaryLambda:list, ImaginaryPartBySqrtm, RealPart, MU); local( j, k, numberOfMu:list, localDeleteFlag:list, nfFraction, limitForLambda_j, Capacity); Capacity = mufixsquare^2+100; preliminaryLambda = listcreate(Capacity); localDeleteFlag = listcreate (Capacity); numberOfMu = listcreate (Capacity); j=1; k=0; numberOfSpheres = 0; preliminaryMu = createPreliminaryMu3mod4( mufixsquare); while (j <= length(preliminaryMu) , MU = preliminaryMu[j]; limitForLambda_j = round(sqrt(norm(nfbasistoalg(K, MU +nfeltmul(K,w,MU))))); /* Create Lambda_j := [a,b]~, which will be stored in the list preliminaryLambda. */ for ( a = -2*limitForLambda_j, 2*limitForLambda_j, for ( b = -limitForLambda_j, limitForLambda_j, /* Check if Lambda/Mu =: nfFraction is in the base rectangle: */ nfFraction = nfeltdiv(K, [a,b]~, MU); RealPart = component( nfFraction, 1) -1/2*component( nfFraction, 2); ImaginaryPartBySqrtm = 1/2*component( nfFraction, 2); if( -1/2 <= RealPart, if( RealPart <= 1/2, if( 0 <= ImaginaryPartBySqrtm, if( ImaginaryPartBySqrtm <= 1/2, if ( idealadd(K, MU, [a,b]~) == [1,0;0,1], k = k+1; numberOfSpheres = numberOfSpheres +1; listput( preliminaryLambda, [a,b]~, numberOfSpheres); listput(numberOfMu,j,numberOfSpheres); listput(localDeleteFlag,0,numberOfSpheres); ); ); ); ); ); ); ); j++; k=0; ); if (length(Lambda) == 0, for (a= 1, length(preliminaryLambda), listput(Lambda,preliminaryLambda[a]); listput(Mu,preliminaryMu[numberOfMu[a]]); listput(deleteFlag,0); ); numberOfSpheres=length(Lambda); removeEverywhereBelowOnes(); cleanseListOfSpheres(); ); localDeleteFlag = deleteSpheresBelow( preliminaryMu,preliminaryLambda,numberOfMu,localDeleteFlag); /* write the global lists Lambda, Mu */ for (a= 1, length(preliminaryLambda), if (localDeleteFlag[a] == 0, listput(Lambda,preliminaryLambda[a]); listput(Mu,preliminaryMu[numberOfMu[a]]); listput(deleteFlag,0); ); ); numberOfSpheres=length(Lambda); print("kept ",numberOfSpheres," hemispheres."); listkill(preliminaryMu); listkill(preliminaryLambda); listkill(numberOfMu); listkill(localDeleteFlag); precomputeSphereData(); }; {addhelp(createSpheres3mod4,"Create the list Lambda with Lambda/Mu hemisphere center.");} /* end of the procedure createSpheres3mod4. */ removeEverywhereBelowOnes() = { local( sphere_center, invertedNormOfMu, centerDistance, k); invertedNormOfMu = vector( length(Mu)); sphere_center = vector( length(Mu)); for( j = 1, length(Mu), invertedNormOfMu[j] = sqrt(1/norm(nfbasistoalg(K, Mu[j]))); sphere_center[j] = nfbasistoalg(K, nfeltdiv(K, Lambda[j], Mu[j])); ); for ( j = 1, length(Lambda), k = j+1; while( deleteFlag[j] != 1 && k <= length(Lambda), if ( deleteFlag[k] != 1, centerDistance = sqrt(norm(sphere_center[j] -sphere_center[k])); if ( invertedNormOfMu[k] <= invertedNormOfMu[j], if ( invertedNormOfMu[k] + centerDistance <= invertedNormOfMu[j], /* Hemisphere number j covers hemisphere number k everywhere. */ deleteFlag[k] = 1; ); ,/* else */ if ( invertedNormOfMu[j] + centerDistance <= invertedNormOfMu[k], /* Hemisphere number k covers hemisphere number j everywhere. */ deleteFlag[j] = 1; ); ); ); k++; ); ); }; /* end of the procedure removeEverywhereBelowOnes. */ createLines() = { local(linesOfSphereList:list); local(emptyListForSphere:list, emptyList:list); numberOfLines= 0; emptyList = listcreate(maxNeighboursNumber+2); linesOfSphere = vector(numberOfSpheres); pointsOfSphere = vector(numberOfSpheres); pointsOfLine = vector(numberOfSpheres); /* Run through the pairs of neighbouring hemispheres, */ /* and compute the line projected from their intersection arc. */ for (j = 1, numberOfSpheres, if ( deleteFlag[j] == 0, if ( radiusSquare[j] == 1, /* for the unit hemispheres */ /* allow two extra lines (the spurious ones) */ linesOfSphereList = listcreate(length(neighbours[j]) +3); pointsOfLine[j] = vector(length(neighbours[j])+4); , /* compared to the smaller hemispheres which have at most one spurious line: */ linesOfSphereList = listcreate(length(neighbours[j])+2); pointsOfLine[j] = vector(length(neighbours[j])+3); ); for ( k = 1, length(neighbours[j]), if ( deleteFlag[neighbours[j][k]] == 0, listput(linesOfSphereList, computeLine(j,k)); ); pointsOfLine[j][k] = emptyList; ); linesOfSphere[j] = linesOfSphereList; listkill(linesOfSphereList); emptyListForSphere = listcreate( length(neighbours[j])*maxNeighboursNumber+4); pointsOfSphere[j] = emptyListForSphere; listkill(emptyListForSphere); ); ); listkill(emptyList); if( mIs3mod4, createVerticalSpuriousLinesOfSpheres3mod4(); createHorizontalSpuriousLinesOfSpheres3mod4(); , /* else */ createVerticalSpuriousLinesOfSpheres(); createHorizontalSpuriousLinesOfSpheres(); ); }; {addhelp(createLines,"Create lists linesOfSphere[j] of projections of intersection arcs of hemispheres. They are lines in the complex plane.");} /* end of the procedure createLines. */ createVerticalSpuriousLinesOfSpheres3mod4() = { /* Create manually elements in linesOfSphere, in order to obtain the spurious edges. */ local( RealPart, radius); for (j = 1, numberOfSpheres, if( deleteFlag[j] == 0, RealPart = component(sphereCenter[j],1) -1/2*component(sphereCenter[j],2); radius = sqrt( radiusSquare[j]); /* if the imaginary axis translated by -1/2 passes through the hemisphere number j, */ if( RealPart -radius <= -1/2, /* Put a line through [-1/2, 0]~ and [0,1]~ into linesOfSphere[j]. */ insertLine( [ [-1/2,0]~, [1/2, 1]~], j); ); /* if the imaginary axis translated by +1/2 passes through the hemisphere number j, */ if( RealPart +radius >= 1/2, /* Put a line through [1/2, 0]~ and [1,1]~ into linesOfSphere[j]. */ insertLine( [ [1/2,0]~, [1/2,1]~], j); ); ); ); }; /* end of the procedure createVerticalSpuriousLinesOfSpheres3mod4(). */ createVerticalSpuriousLinesOfSpheres() = { /* Create manually elements in linesOfSphere, in order to obtain the spurious edges. */ local( RealPart, radius); for (j = 1, numberOfSpheres, if( deleteFlag[j] == 0, RealPart = component(sphereCenter[j],1); radius = sqrt( radiusSquare[j]); /* if the imaginary axis passes through the hemisphere number j, */ if( RealPart -radius < 0, /* Put a line through [0, 0]~ and [0, 1]~ into linesOfSphere[j]. */ insertLine( [ [0,0]~, [0, 1]~], j); ); /* if the imaginary axis translated by +1 passes through the hemisphere number j, */ if( RealPart +radius > 1, /* Put a line through [1, 1]~ and [1,1/2]~ into linesOfSphere[j]. */ insertLine( [ [1,1]~, [0,-1/2]~], j); ); ); ); }; /* end of the procedure createVerticalSpuriousLinesOfSpheres(). */ insertLine( Line, indexOfUnitSphere) = { local( auxiliaryList:list, lineNumber); /* Put Line into linesOfSphere[indexOfUnitSphere]. */ auxiliaryList = linesOfSphere[indexOfUnitSphere]:list; listput( auxiliaryList, Line); linesOfSphere[indexOfUnitSphere] = auxiliaryList; listkill( auxiliaryList); /* Insert an empty pointsOfLine list. */ lineNumber = length( linesOfSphere[indexOfUnitSphere]); auxiliaryList = listcreate(maxNeighboursNumber+3); pointsOfLine[indexOfUnitSphere][lineNumber] = auxiliaryList; listkill( auxiliaryList); }; /* end of the procedure insertLine. */ createHorizontalSpuriousLinesOfSpheres3mod4() = { /* Create manually elements in linesOfSphere, in order to obtain the spurious edges. */ local( indexOfUnitSphere); /* Find the unit hemisphere with center [0,0]~. */ indexOfUnitSphere = findUnitSphere([0,0]~); /* Put a line through [0, 0]~ and [1/2, 0]~ into linesOfSphere[indexOfUnitSphere]. */ insertLine( [ [0,0]~, [1/2, 0]~], indexOfUnitSphere); /* Formerly, the point [[0, 0]~,1] was now put manually into */ /* pointsOfLine[indexOfUnitSphere][lineNumber]. */ /* Find the unit hemisphere with center [1,1]~. */ indexOfUnitSphere = findUnitSphere([1,1]~); /* Put a line through [1, 1]~ and [1/2, 1]~ into linesOfSphere[indexOfUnitSphere]. */ insertLine( [ [1,1]~, [-1/2,0]~], indexOfUnitSphere); /* Formerly, the point [[1, 1]~,1] was now put manually into */ /* Find the unit hemisphere with center [0,1]~. */ indexOfUnitSphere = findUnitSphere([0,1]~); /* Put a line through [0, 1]~ and [1/2, 1]~ into linesOfSphere[indexOfUnitSphere]. */ insertLine( [ [0,1]~, [1/2,0]~], indexOfUnitSphere); /* Formerly, the point [[0, 1]~,1] was now put manually into */ /* pointsOfLine[indexOfUnitSphere][lineNumber]. */ numberOfLines = numberOfLines +3; }; /* end of the procedure createHorizontalSpuriousLinesOfSpheres3mod4(). */ createHorizontalSpuriousLinesOfSpheres() = { /* Create manually elements in linesOfSphere, in order to obtain the spurious edges. */ local( indexOfUnitSphere); /* Find the unit hemisphere with center [0,0]~. */ indexOfUnitSphere = findUnitSphere([0,0]~); /* Put a line through [0, 0]~ and [1/2, 0]~ into linesOfSphere[indexOfUnitSphere]. */ insertLine( [ [0,0]~, [1/2, 0]~], indexOfUnitSphere); /* Formerly, the point [[0, 0]~,1] was now put manually into */ /* pointsOfLine[indexOfUnitSphere][lineNumber]. */ /* Find the unit hemisphere with center [0,1]~. */ indexOfUnitSphere = findUnitSphere([0,1]~); /* Put a line through [0, 1]~ and [1/2, 1]~ into linesOfSphere[indexOfUnitSphere]. */ insertLine( [ [0,1]~, [1/2, 0]~], indexOfUnitSphere); /* Formerly, the point [[0, 1]~,1] was now put manually into */ /* pointsOfLine[indexOfUnitSphere][lineNumber]. */ /* Find the unit hemisphere with center [1,1]~. */ indexOfUnitSphere = findUnitSphere([1,1]~); /* Put a line through [1, 1]~ and [1/2, 1]~ into linesOfSphere[indexOfUnitSphere]. */ insertLine( [ [1,1]~, [-1/2,0]~], indexOfUnitSphere); /* Formerly, the point [[1, 1]~,1] was now put manually into */ /* pointsOfLine[indexOfUnitSphere][lineNumber]. */ /* Find the unit hemisphere with center [1,0]~. */ indexOfUnitSphere = findUnitSphere([1,0]~); /* Put a line through [1, 0]~ and [1/2, 0]~ into linesOfSphere[indexOfUnitSphere]. */ insertLine( [ [1,0]~, [-1/2,0]~], indexOfUnitSphere); /* Formerly, the point [[1, 0]~,1] was now put manually into */ /* pointsOfLine[indexOfUnitSphere][lineNumber]. */ numberOfLines = numberOfLines +4; }; /* end of the procedure createHorizontalSpuriousLinesOfSpheres(). */ removePoints() = { /* Wipe points which are already deleted in the pointsOfSphere lists out of the pointsOfLine lists. */ local(newPointsList:list, point); for (r = 1, numberOfSpheres, for ( s = 1, length(pointsOfLine[r]), newPointsList = listcreate(1+length(pointsOfLine[r][s]:list)); for ( t = 1, length(pointsOfLine[r][s]:list), point = pointsOfLine[r][s]:list[t]; for (j = 1, numberOfSpheres, if (setsearch(Set(pointsOfSphere[j]:list), point), listput(newPointsList, point); listsort(newPointsList,1); ); ); ); pointsOfLine[r][s] = newPointsList; listkill(newPointsList); ); ); }; {addhelp(removePoints,"Wipe points which are already deleted in the pointsOfSphere lists out of the pointsOfLine lists.");} /* end of procedure removePoints. */ removeLines() = { /* Remove lines which meet the surface at less than two vertices. */ local(newLinesList: list, newPointsOfLineList: list); for (r = 1, numberOfSpheres, newLinesList = listcreate( length( linesOfSphere[r]:list)); newPointsOfLineList = listcreate( length( pointsOfLine[r])); for ( s = 1, length( pointsOfLine[r]), /* Keep only lines with at least two visible points: */ if (length(pointsOfLine[r][s]) >= 2, listput( newLinesList, linesOfSphere[r]:list[s]); listput( newPointsOfLineList, pointsOfLine[r][s]); ); ); linesOfSphere[r] = newLinesList; pointsOfLine[r] = Vec(newPointsOfLineList); listkill(newLinesList); listkill(newPointsOfLineList); ); }; {addhelp(removeLines,"Remove lines which meet the surface at less than two vertices.")}; /* end of the procedure removeLines */ getVertices() = { local( pointAtIntersection, RealPart, ImaginaryPartBySqrtm); for ( j:small = 1, numberOfSpheres, for ( k:small = 1, length(linesOfSphere[j]:list), for ( r:small = 1, length(linesOfSphere[j]:list), if(linesOfSphere[j]:list[k] != linesOfSphere[j]:list[r], pointAtIntersection = intersectLines(j,k,r); if( mIs3mod4, 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, addVertex( pointAtIntersection, j); addVertexToLine( pointAtIntersection, j, k); addVertexToLine( pointAtIntersection, j, r); ); ); ); ); ,/* else m not congruent to 3 modulo 4 */ if( /* pointAtIntersection is in the base rectangle */ 0 <= component(pointAtIntersection,1) && component(pointAtIntersection,1) <= 1 && 0 <= component(pointAtIntersection,2) && component(pointAtIntersection,2) <= 1, addVertex( pointAtIntersection, j); addVertexToLine( pointAtIntersection, j, k); addVertexToLine( pointAtIntersection, j, r); ); ); ); ); ); ); /* Get 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)); ); /*Could be done now: Delete points which are below another point in vertexList */ }; /* end of procedure getVertices */ cleanseVertexLists( deletePoint) = { local( newVertexList: list); /* Clean the set of vertices (totalPointSet) */ newVertexList = listcreate(length(totalPointSet)); for ( j = 1, length(totalPointSet), if ( deletePoint[j] == 0, listput(newVertexList, eval(totalPointSet[j])); ); ); totalPointSet = Set( newVertexList); /* Clean the lists pointsOfSphere[k] of vertices by throwing out the ones that are no more in totalPointSet */ for ( k = 1, numberOfSpheres, listkill( newVertexList); for ( l = 1, length( pointsOfSphere[k]:list), /* Check that pointsOfSphere[k][l] still exists in totalPointSet */ if ( setsearch( totalPointSet, pointsOfSphere[k]:list[l]), listput( newVertexList, pointsOfSphere[k]:list[l]); ); ); pointsOfSphere[k] = newVertexList; ); listkill( newVertexList); }; /* end of procedure cleanseVertexLists */ compareHeights( numberOfRequiredCorners) = { local(z, comparisonHeightSquare, heightSquare, setOfSingularPoints); local(deletePoint, numberOfErased, point, auxiliaryList:list); numberOfErased = 0; deletePoint = vector( length(totalPointSet)); listsort( singularpoint:list); setOfSingularPoints = Set(singularpoint:list); /* z runs through the projections of intersection points */ for ( j = 1, length(totalPointSet), point = eval(totalPointSet[j]); z = component( point, 1); if( /* z not a singular point */ setsearch(setOfSingularPoints,z,1), /* then find the hemispheres touching the surface above z */ comparisonHeightSquare = component( point, 2); if (comparisonHeightSquare <= 0, deletePoint[j] = 1; , /* else consider z */ for (k = 1, numberOfSpheres, if(deleteFlag[k] == 0, /* Height square of hemisphere k over the point z: */ heightSquare = radiusSquare[k] -norm( nfbasistoalg( K, z -sphereCenter[k] )); if( heightSquare > comparisonHeightSquare, /* Hemisphere number k lies above the point z. */ /* delete z */ deletePoint[j] = 1; ); ); ); ); ); ); /* Throw out the vertices which are marked for erasure by the vector deletePoint. */ cleanseVertexLists( deletePoint); if( mIs3mod4, /* Now put the remaining vertices into all of the lists of the hemispheres on which they lie. */ /* z runs through the projections of intersection points */ for ( j = 1, length(totalPointSet), point = eval(totalPointSet[j]); z = component( point, 1); comparisonHeightSquare = component( point, 2); /* Find the hemispheres touching z: */ for (k = 1, numberOfSpheres, if(deleteFlag[k] == 0, /* Height square of hemisphere k over the point z: */ heightSquare = radiusSquare[k] -norm(nfbasistoalg( K, z -sphereCenter[k] )); if( heightSquare == comparisonHeightSquare, /*Check if (z,heightSquare) is already in pointsOfSphere[k] */ if( setsearch( Set(pointsOfSphere[k]:list), totalPointSet[j]), ,/* Else put (z,heightSquare) into pointsOfSphere[k]. */ print("Located ",totalPointSet[j]," on the hemisphere centered at ", sphereCenter[k]); auxiliaryList = pointsOfSphere[k]:list; listput(auxiliaryList, point); pointsOfSphere[k]:list = auxiliaryList; listkill(auxiliaryList); ); ); ); ); ); ); if ( numberOfRequiredCorners > 0, for( k = 1, numberOfSpheres, if( length( pointsOfSphere[k]:list) < numberOfRequiredCorners, deleteFlag[k] = 1; /* Delete the hemisphere with center Lambda[k]/Mu[k]. */ numberOfErased++; ); ); print( numberOfErased, " hemispheres were deleted because they touched the surface less than ", numberOfRequiredCorners," times."); ); cleanseListOfSpheres(1); /* option keepPoints selected. */ }; /* end of procedure compareHeights. */ getLinesAndVertices() = { logfile = Str( "m=",m," Relevant hemispheres in hyperbolic space found in ",gettime()," milliseconds."); print("Computing intersection lines of ",numberOfSpheres, " hemispheres..."); createLines(); print("... recorded ",numberOfLines," lines, ..."); getVertices(); print("... which have ",length(totalPointSet)," vertices."); logfile = Str( logfile,"m=",m," Lines and vertices computation took ",gettime()," milliseconds."); compareHeights(0); logfile = Str( logfile,"m=",m,"Killing low vertices took ", gettime()," milliseconds."); }; /* end of procedure getLinesAndVertices */