\\ --------------- 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 1.0.1 of October 31st, 2010. \\--------------------------------------------------------------------------------------- conjugate(y) = nfalgtobasis(K,conj(nfbasistoalg(K,y))); PoincareAction(a,b,c,d,z,zetasquare) = { /* Compute Poincare's action of the matrix [a,b;c,d] on the point (z,zeta) */ /* in Hyperbolic Space, using Floege's version of Poincare's formula. */ local(dMINUScz, divised, divisor, cConjugated, dMINUSczConjugated); if( type(b) != "t_COL", /* if the matrix was entered by hand, adjust type of the entries. */ b = nfalgtobasis(K,b); c = nfalgtobasis(K,c); d = nfalgtobasis(K,d); if( type(z) != "t_COL", z = nfalgtobasis(K,z); ); ); dMINUScz = -nfeltmul(K,c,z) +d; dMINUSczConjugated = conjugate(dMINUScz); cConjugated = conjugate(c); divised = nfeltmul(K, nfeltmul(K,a,z)-b, dMINUSczConjugated ) -zetasquare*nfeltmul(K,a,cConjugated); divisor = nfeltmul(K, dMINUScz, dMINUSczConjugated) +zetasquare*nfeltmul(K,c,cConjugated); if( component( divisor, 2) == 0, if( type(a) != "t_COL", /* if the matrix was entered by hand */ print("height square ",zetasquare/component(divisor,1)^2," over"); ); , /* else the computed height will not be a real number: contradiction. */ print("Contradiction: ",divisor," not real, so the height will not be a real number."); ); /* Return */ nfeltdiv(K, divised, divisor) }; addhelp(PoincareAction, "Compute the translate under the matrix [a, b; c, d], with respect to the number field K, of a point (z, Zeta) in Hyperbolic space. Function called as PoincareAction( a,b,c,d,z,zetasquare). Use K = nfinit(w^2+m) before, if necessary."); /* end of the function PoincareAction. */ recordOrbitSum() = { local(orbitSum); orbitSum = vector(numberOfSpheres +m); for (j = 1, numberOf2cells, for( r = 1, length( cornersOf2cell[j]), /* Write down a binary information about which vertex orbits are represented among the elements of cornersOf2cell[j]. */ orbitSum[j] = orbitSum[j] + 2^vertexOrbitNumber[cornersOf2cell[j][r]]; ); ); /* return */ orbitSum }; /* end of function orbitSum */ computeSingularStabilizer(s)= { local(a,b,c,d, PointAndDenominatorCouple); local(currentStabilizer:list,numberInSet, c_max, potential_c, cs, csSquare); currentStabilizer = listcreate(44); /* With c = 0, we get just the trivial stabilizing matrices (except for m= 1 or 3): */ listput(currentStabilizer, matid(2) ); listput(currentStabilizer, -matid(2) ); if (m == 1, c = [0,0]~; for ( k = 0, 1, d = [0, (-1)^k]~; a = nfalgtobasis(K, conj( nfbasistoalg(K,d))); b = nfeltmul(K, a-d, s); listput(currentStabilizer, [a,b; c,d] ); ); ); /* Now we search stabilizer elements with nonzero entry c. */ numberInSet = setsearch(Set(singularpoint), s); PointAndDenominatorCouple = eval(Set(singularDenominator)[numberInSet]); if ( s != component(PointAndDenominatorCouple, 1), error("Set function gives wrong order in the singular points list. Check the function computeSingularStabilizer"); ); c_max = component(PointAndDenominatorCouple, 2); /* Use the function createPreliminaryMu to get a list of all integers in the number field, */ /* which have norm at most c_max. */ potential_c = Set(); for( j = 1, c_max^2, if( mIs3mod4, potential_c = setunion( potential_c, Set(createPreliminaryMu3mod4(j))); , /* else */ potential_c = setunion( potential_c, Set(createPreliminaryMu(j))); ); ); while( length( currentStabilizer) < 23, /* We record only 23 of the infinitely many elements in the stabilizer. */ for (j = 1, length(potential_c), c = eval(potential_c[j]); /* Check the two conditions that c is such that cs and cs^2 are integers */ /* in the number field. This is an ideal in the ring of integers. */ cs = nfeltmul( K, c, s); csSquare = nfeltmul(K, cs, s); if( frac(cs) == [0,0]~ && frac(csSquare) == [0,0]~, /* print("singular point ",s, */ /* " with stabilizer the free abelian group ", */ /* "generated by [",c,"]_1 = " ); */ a = [1,0]~ -cs; b = -csSquare; d = [1,0]~ +cs; if( length( currentStabilizer) < 23, listput(currentStabilizer, [a,b; c,d] ); ); if( length( currentStabilizer) < 23, /* Record the inverse [d,-b;-c,a] of [a,b;c,d]. */ listput(currentStabilizer, [d,-b;-c,a] ); ); /* "record also [",c,"]_-1 = " */ a = [-1,0]~ -cs; d = [-1,0]~ +cs; if( length( currentStabilizer) < 23, listput(currentStabilizer, [a,b; c,d] ); ); if( length( currentStabilizer) < 23, /* Record the inverse [d,-b;-c,a] of [a,b;c,d]. */ listput(currentStabilizer, [d,-b;-c,a] ); ); listsort( currentStabilizer, 1); /* Flag 1 to delete double occurencies. */ ); ); c_max++; for( j = (c_max -1)^2 +1, c_max^2, if( mIs3mod4, potential_c = setunion( potential_c, Set(createPreliminaryMu3mod4(j))); , /* else */ potential_c = setunion( potential_c, Set(createPreliminaryMu(j))); ); ); ); /* Return */ currentStabilizer }; /* end of the function computeSingularStabilizer. */ computeVertexStabilizer(z, rsquare)= { local(a,b,c,d, q, sBound, cz_real, cz_omegacoeff, discriminant, squareRoot); local(currentStabilizer:list, kBound, jBound, rsquareInverse, cz, cNormsquare, rcNormsquare); currentStabilizer = listcreate(24); /* The maximal order of finite subgroups in the Bianchi groups is 24. */ rsquareInverse = 1/rsquare; /* With c = 0, we get just the trivial stabilizing matrices (except for m= 1 or 3): */ listput(currentStabilizer, matid(2) ); listput(currentStabilizer, -matid(2) ); if (m == 1, c = [0,0]~; for ( k = 0, 1, d = [0, (-1)^k]~; a = nfalgtobasis(K, conj( nfbasistoalg(K,d))); b = nfeltmul(K, a-d, z); listput(currentStabilizer, [a,b; c,d] ); ); ); /* Now we search stabilizer elements with nonzero entry c. */ /* c = j +k*w runs through the lattice of integers with 0 < |c| <= 1/r. */ /* The equation |c| = sqrt(j^2 +k^2*m) gives us the following bound for k. */ kBound = 1/(sqrt(rsquare*m)); for (k = ceil(-kBound), floor(kBound), jBound = sqrt(rsquareInverse -m*k^2); for (j = ceil(-jBound), floor(jBound), c = nfalgtobasis(K,j + k*w); cNormsquare = norm(nfbasistoalg(K,c)); if ( cNormsquare > rsquareInverse, print("***Error in function computeVertexStabilizer on stabilizer entry c."); ); if ( c != [0,0]~ , /* d runs through the lattice of integers with |cz-d|^2 +r^2*|c|^2 = 1. Use d =: q +s*w. */ cz = nfeltmul( K, c, z); cz_real = component( cz, 1); /* check that 'a' is in the ring of integers */ if( frac( 2*cz_real) == 0, cz_omegacoeff = component( cz, 2); rcNormsquare = rsquare*cNormsquare; sBound = abs(cz_omegacoeff) +sqrt((1 -rcNormsquare)/m); for ( s = -ceil(sBound), ceil(sBound), discriminant = 1 -rcNormsquare -m*(cz_omegacoeff -s)^2; if( issquare( discriminant, &squareRoot), /* q+ */ q = round( cz_real + squareRoot); if( frac(q) == 0, d = nfalgtobasis(K, q+ s*w); if (norm(nfbasistoalg(K, cz -d)) +rcNormsquare == 1 && idealadd(K,c,d) == 1, a = nfalgtobasis(K, conj( nfbasistoalg(K,d))) -2*cz_real*[1,0]~; b = nfeltdiv(K, nfeltmul(K,a,d) -[1,0]~, c); /* Check that b is in the ring of integers: */ if ( frac(b) == [0,0]~, listput(currentStabilizer, [a,b; c,d] ); ); ); ); ); if( discriminant > 0, /* q- */ q = round( cz_real - squareRoot); if( frac(q) == 0, d = nfalgtobasis(K, q+ s*w); if ( norm(nfbasistoalg(K, cz -d)) +rcNormsquare == 1 && idealadd(K,c,d) == 1, a = nfalgtobasis(K, conj( nfbasistoalg(K,d))) -2*cz_real*[1,0]~; b = nfeltdiv(K, nfeltmul(K,a,d) -[1,0]~, c); /* Check that b is in the ring of integers: */ if ( frac(b) == [0,0]~, listput(currentStabilizer, [a,b; c,d] ); ); ); ); ); ); ); ); ); ); /* Return */ currentStabilizer }; /* end of function computeVertexStabilizer */ computeVertexStabilizer3mod4( z, rsquare)= { local(a,b,c,d, s_limitSummand, cz, cz_R, cz_W, discriminant, squareRoot); local( cNorm, k_limitSummand, j_limit, currentStabilizer:list ); currentStabilizer = listcreate(41); /* With c = 0, we get just the trivial stabilizing matrices (except for m= 1 or 3): */ listput(currentStabilizer, matid(2) ); listput(currentStabilizer, -matid(2) ); /* Now we search stabilizer elements with nonzero entry c. */ /* c = j +k*w runs through the lattice of integers with 0 < |c| <= 1/r */ j_limit = floor(sqrt((1+1/m)/rsquare)); /* j runs from -j_limit through j_limit: */ for (j = -j_limit, j_limit, /* |c| = sqrt( (j -k/2)^2 +m*(k/2)^2) = sqrt( j^2 +(m+1)/4*k^2 -j*k) */ /* therefore, k_limit(+,-) = 2*j/(m+1) (+,-)2*sqrt(-j^2*m +(m+1)/(rsquare))/(m+1) */ discriminant = -j^2*m +(m+1)/(rsquare); if( discriminant >= 0, k_limitSummand = 2*sqrt(discriminant)/(m+1); for (k = floor( 2*j/(m+1) - k_limitSummand), ceil( 2*j/(m+1) + k_limitSummand), c = nfalgtobasis(K,j + k*w); /* Check that |c| <= 1/r and c nonzero. */ if ( norm(nfbasistoalg(K,c)) <= 1/rsquare && c != [0,0]~ , /* d runs through the lattice of integers with |cz-d|^2 +r^2|c|^2 = 1. */ /* Decompose d as q +s*w. */ cz = nfeltmul( K, c, z); cz_R = component( cz, 1); cz_W = component( cz, 2); cNorm = norm(nfbasistoalg(K,c)); /* s_limit(+,-) = cz_W (+,-) 2*sqrt( (1 -r^2|c|^2)/m). */ s_limitSummand = 2*sqrt( (1 -rsquare*cNorm)/m); for ( s = floor(cz_W -s_limitSummand), ceil(cz_W +s_limitSummand), discriminant = 1 -rsquare*cNorm -(m/4)*cz_W^2 +(m/2)*cz_W*s -(m/4)*s^2; if( issquare( discriminant, &squareRoot), /* q(+,-) = cz_R -cz_W/2 +s/2 (+,-) sqrt( discriminant) */ qPlus = cz_R -cz_W/2 +s/2 +squareRoot; if( frac( qPlus) == 0, currentStabilizer = concat(currentStabilizer, FinishStabilizer( round(qPlus),s,c,cNorm,cz, cz_R, cz_W, z,rsquare)); ); if( discriminant != 0, qMinus = cz_R -cz_W/2 +s/2 -squareRoot; if( frac( qMinus) == 0, currentStabilizer = concat(currentStabilizer, FinishStabilizer( round(qMinus),s,c,cNorm,cz, cz_R, cz_W, z,rsquare)); ); ); ); ); ); ); ); ); /* Return */ currentStabilizer }; /* end of function computeVertexStabilizer3mod4. */ FinishStabilizer(q,s,c,cNorm,cz, cz_R, cz_W, z,rsquare) = { /* Finish the computation of the stabilizer of (z,rsquare) in Hyperbolic Space. */ /* m congruent 3 modulo 4. */ local(a,b,d,StabilizingMatrix:list, ad, TwoRe_cz); StabilizingMatrix = listcreate(1); /* Floege deduces a = conj(d) -2*Re(cz). */ TwoRe_cz = 2*cz_R -cz_W; /* Check that a is in the ring of integers: */ if( frac(TwoRe_cz) == 0, d = nfalgtobasis(K, q + s*w); /* Check that |cz-d|^2 +r^2*|c|^2 == 1, */ if( norm(nfbasistoalg(K, cz -d)) +rsquare*cNorm == 1 /* and check that c,d generate the ring of integers. */ && idealadd(K,c,d) == 1, a = conjugate(d) -TwoRe_cz*[1,0]~; ad = nfeltmul(K,a,d); /* b is determined by the determinant 1 equation: */ b = nfeltdiv(K, ad -[1,0]~, c); /* Check that b is in the ring of integers: */ if ( frac(b) == [0,0]~, /* Check that z == conj(d-cz)(az-b) -r^2*conj(c)a. */ if( z == nfeltmul(K, conjugate(d -cz), nfeltmul(K,a,z) -b) -rsquare*nfeltmul(K, conjugate(c), a), /* Check stabilization. */ if( z != PoincareAction(a,b,c,d,z,rsquare), print("***Error in function FinishStabilizer ", "for the stabilization of z = ",z, ", rsquare = ",rsquare,"\n [a,b; c,d](z,rsquare) = ", PoincareAction(a,b,c,d,z,rsquare) ); ); listput(StabilizingMatrix, [a,b; c,d] ); , /* else */ error(Str("***Error in function FinishStabilizer ", "caused by [a,b; c,d] = ",[a,b; c,d]) ); ); ); ); ); /* Return */ StabilizingMatrix }; /* end of subfunction FinishStabilizer of the function computeVertexStabilizer3mod4. */ getIdentificationMatrices(p, p_2, stabilizerCardinal) = { /* For m is not 1, and not congruent to 3 modulo 4: */ /* Computation of the identification matrices [a,b; c,d] */ /* which transport the point p=(z,r) into the point p_2 =(Zeta, rho)*/ /* in hyperbolic space, where r and rho are both positive heights.*/ /* If p and p_2 are cusps, use the function getSingularIdentifications later. */ local(z,rsquare, Zeta, rhosquare,a,b,c,d, hBound, cz, cz_real, cz_omegacoeff, kbound, dConj, rByRho, oneByRrho, squareRoot); local(vertexTransport:list, g, discriminant, jbound, lefthandside, rsquareByRhosquare, rZetaByRhoSquare); vertexTransport = listcreate(42); z = component(p,1); rsquare = component(p,2); Zeta = component(p_2,1); rhosquare = component(p_2,2); if ( rsquare > 0 && rhosquare > 0, rsquareByRhosquare = rsquare/rhosquare; if( issquare(rsquareByRhosquare, &rByRho), ,/* else */ rByRho = sqrt(rsquareByRhosquare) ); rZetaByRhoSquare = rsquareByRhosquare*nfeltmul(K,Zeta,Zeta); /* First case: c = 0 */ c = [0,0]~; for( n = 0, 1, d = [(-1)^n, 0]~; a = d; /* = nfeltdiv(K,1,d) */ b = nfeltmul(K, a, z) - nfeltmul(K, d, Zeta); dConj = conj(nfbasistoalg(K,d)); /* check that b is in the ring of integers: */ if ( frac(b) == [0,0]~, /* check that [1,-z-Zeta; 0,1]p = p_2:*/ lefthandside = nfeltmul(K, dConj, (nfeltmul(K,a,z) -b)); if( /* (norm(nfbasistoalg(K,d)))^2 = */ 1 == rsquareByRhosquare && nfeltmul(K,lefthandside,lefthandside) == rZetaByRhoSquare, /* [a,b; 0,d] sends p to p_2 */ listput(vertexTransport, [a,b; c,d]); ); ); ); /* Second case: c is nonzero.*/ /* c = j +kw runs through the ring of integers of the number field, verifying */ /* 1/(r*rho) >= normsquare(c) = j^2 +mk^2 */ oneByRrho = rByRho/rsquare; kbound = floor( sqrt( oneByRrho/m)); for ( k = -kbound, kbound, jbound = round( sqrt( oneByRrho -m*k^2)); for( j = -jbound, jbound, c = j +k*w; if( c != 0 && (norm(c))^2 <= 1/(rsquare*rhosquare), cz = nfeltmul( K, c, z); cz_real = component( cz, 1); cz_omegacoeff = component( cz, 2); cz = nfbasistoalg(K, cz); /* d =: g+hw runs through the ring of integers of the number field, verifying */ /* normsquare(cz -d) +r^2*normsquare(c) = r/rho. */ hBound = abs(cz_omegacoeff) +sqrt( (rByRho -rsquare*norm(c))/m); /* hBound >= abs(h) */ for( h = -round(hBound), round( hBound), discriminant = rByRho -m*(cz_omegacoeff-h)^2 -rsquare*norm(c); /* Check that 'discriminant' is a rational square: */ if( type(discriminant) == "t_FRAC" || type(discriminant) == "t_INT", if( issquare( discriminant, &squareRoot), /* obtain g with the positive squareroot */ g = cz_real +squareRoot; /* Check that g is a rational integer */ if( frac(g) == 0, vertexTransport = concat( vertexTransport, getRemainingEntries( round(g),h,c,z,cz,Zeta,rsquare,rhosquare, rsquareByRhosquare, rZetaByRhoSquare, rByRho)); ); /* obtain g with the negative squareroot */ g = cz_real -squareRoot; /* Check that g is a rational integer */ if( frac(g) == 0, vertexTransport = concat( vertexTransport, getRemainingEntries( round(g),h,c,z,cz,Zeta,rsquare,rhosquare, rsquareByRhosquare, rZetaByRhoSquare, rByRho)); ); ); ); ); ); ); ); if( length(vertexTransport) > 0, /* Delete double occurencies in the vertexTransport list: */ listsort(vertexTransport, 1); if( abs(z -Zeta) != [1,0]~ && abs(z -Zeta) != [0,1]~, /* There are length(vertexTransport) matrices sending p to p_2. */ ); if( length( vertexTransport) != stabilizerCardinal, print("***Error in function getIdentificationMatrices, on vertexTransport from ", p," to ",p_2,": ","It must be a coset of the two vertex stabilizers!"); for (r = 1, length(vertexTransport), print(vertexTransport[r]); print(); ); ); ); ); /* return the list */ vertexTransport }; /* end of function getIdentificationMatrices */ getRemainingEntries(g,h,c,z,cz,Zeta,rsquare,rhosquare, rsquareByRhosquare, rZetaByRhoSquare, rByRho) = { local(a,b,d, d_minus_czConj, lefthandside, TransportMatrix:list); TransportMatrix = listcreate(1); d = g +h*w; d_minus_czConj = conj(d -cz); /* Check the square of the equation |cz -d|^2 +r^2|c|^2 = r/rho to be verified by d: */ if( (norm(cz -d) +rsquare*norm(c))^2 == rsquareByRhosquare, a = nfalgtobasis( K, 1/rByRho*d_minus_czConj) -nfeltmul(K, c, Zeta); /* Check that the matrix entry "a" has integer entries: */ if( frac(a) == [0,0]~, a = round(a); /* Get b by the determinant 1 equation: */ b = nfeltdiv(K, nfeltmul(K, a, d) -[1,0]~, c); /* check that b has integer entries: */ if( (frac(component(b,1)))^3 == 0 && (frac(component(b,2)))^3 == 0, b = round(b); /* Check that [a,b; c,d]p == p_2 by the following equation: */ /* nfeltmul(K, conj(d -cz), (nfeltmul(K,a,z)-b)) -rsquare*nfeltmul(K,conj(c),a) == r*Zeta/rho */ lefthandside = nfeltmul(K, d_minus_czConj, (nfeltmul(K,a,z) -b)) -rsquare*nfeltmul(K, conj(c), a); if( lefthandside == Zeta*rByRho, c = nfalgtobasis(K, c); d = nfalgtobasis(K, d); /* [a,b; c,d] sends p to p_2 */ listput(TransportMatrix, [a,b; c,d]); ); ); ); ); TransportMatrix }; /* end of the common subfunction getRemainingEntries of the functions getIdentificationMatrices and getIdentificationMatrices3mod4.*/ getIdentificationMatrices3mod4(p, p_2, stabilizerCardinal) = { /* For m is not 3, but congruent to 3 modulo 4: */ /* Computation of the identification matrices [a,b; c,d] */ /* which transport the point p=(z,r) into the point p_2 =(Zeta, rho)*/ /* in hyperbolic space, i.e., r and rho are both positive heights.*/ /* If p and p_2 are cusps, use the function getSingularIdentifications later. */ /* Check that the cardinal of the stabilizer is the same for p and p_2 */ /* before calling this function. */ local(z,rsquare, Zeta, rhosquare,a,b,c,d, hBound, cz, cz_real, cz_omegacoeff, rZetaByRhoSquare, squareRoot); local(vertexTransport:list, g, discriminant, jbound, k_limitPlus, k_limitMinus, rsquareByRhosquare); local(lefthandside, ThereIsNoRamification, s_limitMinus, s_limitPlus,q, rByRho, rRhoInverse, rcSquare); vertexTransport = listcreate(24); z = component(p,1); rsquare = component(p,2); Zeta = component(p_2,1); rhosquare = component(p_2,2); if ( rsquare > 0 && rhosquare > 0, rsquareByRhosquare = rsquare/rhosquare; if( issquare(rsquareByRhosquare, &rByRho), ,/* else */ rByRho = sqrt(rsquareByRhosquare) ); rZetaByRhoSquare = rsquareByRhosquare*nfeltmul(K,Zeta,Zeta); /* First case: c = 0 */ c = [0,0]~; /* d = g+hw = g +h(-1/2 +1/2sqrt(-m)) runs through the ring of integers with |d|^2 = r/rho. */ /* |d|^2 = (g -h/2)^2 + m(h/2)^2 */ /* |d|^2 = g^2 -gh +h^2/4 +mh^2/4 */ hBound = sqrt( 4*rByRho/(m+1)); for( h = -ceil( hBound), ceil( hBound), /* g_+- = h/2 +- sqrt( r/rho -mh^2/4). */ /* If the above discriminant is zero, run just one case, else two cases for g_+-. */ if( rsquareByRhosquare == (m*h^2/4)^2, ThereIsNoRamification = 1; , /* else */ ThereIsNoRamification = 0; ); for ( Case = ThereIsNoRamification, 1, discriminant = rByRho -m*h^2/4; if( issquare( discriminant, &squareRoot), g = h/2 +(-1)^Case*squareRoot; /* Check that g is a rational integer */ if( frac(g) == 0, d = nfalgtobasis(K, round(g)+h*w); /* Check that |d|^2 = r/rho */ if( (norm(nfbasistoalg(K,d)))^2 == rsquareByRhosquare, a = nfeltdiv(K,1,d); /* check that 'a' is in the ring of integers: */ if ( frac(a) == [0,0]~, b = nfeltmul(K, a, z) - nfeltmul(K, d, Zeta); /* check that b is in the ring of integers: */ if ( frac(b) == [0,0]~, /* check that [a,b; 0,d]p = p_2:*/ lefthandside = nfeltmul(K, conj(nfbasistoalg(K,d)), (nfeltmul(K,a,z) -b)); if( nfeltmul(K,lefthandside,lefthandside) == rZetaByRhoSquare, /* [a,b; 0,d] sends p to p_2 */ listput(vertexTransport, [a,b; c,d]); ); ); ); ); ); ); ); ); /* Second case: c is nonzero.*/ /* The equation (2') gives us $ r^2|c|^2 \leq \frac{r}{\rho} $. */ /* c = j +kw runs through the ring of integers of the number field, verifying */ /* 1/(r*rho) >= normsquare(c) > 0. */ rRhoInverse = rByRho/rsquare; jbound = ceil( sqrt( (1+1/m)*rRhoInverse )); /* We find the extremal values of j by an argument analogous to the one in the documentation of computeVertexStabilizer3mod4. */ for( j = -jbound, jbound, /* We have |c|^2 = (j -k/2)^2 +m(k/2)^2 = j^2 +((m+1)k^2)/4 -jk, */ /* therefore k_limit^2 + 4/(m+1)*(-j*k_limit +j^2 -1/(r*rho)) = 0. */ if( (m+1)*rRhoInverse >= j^2*m, /* (check discriminant >= 0 ) */ /* Thus let */ k_limitPlus = 2/(m+1)*j +2*sqrt(( (m+1)*rRhoInverse -j^2*m)/(m+1)); k_limitMinus = 4/(m+1)*j -k_limitPlus; /* = 2/(m+1)*j -2*sqrt(( (m+1)*rRhoInverse -j^2*m)/(m+1)). */ for ( k = floor(k_limitMinus), ceil(k_limitPlus), c = j +k*w; if( c != 0 && (norm(c))^2*(rsquare*rhosquare) <= 1, /* Decompose cz as cz_real +cz_omegacoeff*w with rational integer coefficients. */ cz = nfeltmul( K, c, z); cz_real = component( cz, 1); cz_omegacoeff = component( cz, 2); cz = nfbasistoalg(K, cz); rcSquare = rsquare*norm(c); /* d =: q+s*w runs through the ring of integers of the number field, verifying */ /* normsquare(cz -d) +r^2*normsquare(c) = r/rho. */ /* We decompose d as q +sw, where q and s are rational integers. */ s_limitPlus = cz_omegacoeff + 2*sqrt( (rByRho -rcSquare )/m); s_limitMinus = 2*cz_omegacoeff -s_limitPlus; /* = cz_omegacoeff - 2*sqrt( (rByRho -rsquare*norm(c))/m) */ for( s = floor(s_limitMinus), ceil(s_limitPlus), discriminant = rByRho -rcSquare -m/4*(cz_omegacoeff -s)^2 ; /* Check that 'discriminant' is a rational square: */ if( type(discriminant) == "t_FRAC" || type(discriminant) == "t_INT", if( issquare( discriminant, &squareRoot), /* obtain q+ with the positive squareroot */ q = cz_real -cz_omegacoeff/2 +s/2 +squareRoot; if( frac(q) == 0, /* Check that q is a rational integer */ vertexTransport = concat( vertexTransport, getRemainingEntries(round(q),s,c,z,cz,Zeta,rsquare,rhosquare, rsquareByRhosquare, rZetaByRhoSquare, rByRho) ); ); /* obtain q- with the negative squareroot */ q = cz_real -cz_omegacoeff/2 +s/2 -squareRoot; if( frac(q) == 0, /* Check that q is a rational integer */ vertexTransport = concat( vertexTransport, getRemainingEntries( round(q),s,c,z,cz,Zeta,rsquare,rhosquare, rsquareByRhosquare, rZetaByRhoSquare, rByRho) ); ); ); ); ); ); ); ); ); if( length(vertexTransport) > 0, /* Delete double occurencies in the vertexTransport list: */ listsort(vertexTransport, 1); if( abs(z -Zeta) != [1,0]~ && abs(z -Zeta) != [0,1]~, /* There are length(vertexTransport) matrices sending p to p_2. */ ); if( length( vertexTransport) != stabilizerCardinal, error(Str("***Error in function getIdentificationMatrices3mod4, on vertexTransport from ", p," to ",p_2,": ","It must be a coset of the two vertex stabilizers!")); for (r = 1, length(vertexTransport), print(vertexTransport[r]); print(); ); ); ); ); /* return the list */ vertexTransport }; /* end of function getIdentificationMatrices3mod4 */ recordVectorOfIdentifications() = { local( vectorOfKnownMatrices, setOfKnownMatrices, vectorOfVectors); /* Write all matrices which occured as identification matrices into a vector. */ vectorOfVectors = Vec( vertexIdentifications); vectorOfKnownMatrices = []; for( k = 1, length(vectorOfVectors), vectorOfKnownMatrices = concat(vectorOfKnownMatrices, Vec(vectorOfVectors[k])); ); setOfKnownMatrices = Set(); for( k = 1, length(vectorOfKnownMatrices), setOfKnownMatrices = setunion(setOfKnownMatrices, vectorOfKnownMatrices[k]); ); vectorOfKnownMatrices = vector( length( setOfKnownMatrices)); for( k = 1, length( setOfKnownMatrices), vectorOfKnownMatrices[k] = eval(setOfKnownMatrices[k]); ); /* return */ vectorOfKnownMatrices; }; /* end of function recordVectorOfIdentifications */ checkIdentification(origin_r, end_r, origin_s, end_s) = { /* We intersect the set of matrices transporting */ /* the origin of the edge r to the origin of the edge s */ /* with the one transporting the end of the edge r to the end of the edge s. */ setintersect( vertexIdentifications[origin_r, origin_s], vertexIdentifications[end_r, end_s] ) /* The matrices in this intersection set are sending eval(totalPointSet[origin_r]) to eval(totalPointSet[origin_s]) and eval(totalPointSet[end_r]) to eval(totalPointSet[end_s]) forth and back. */ /* Return this intersection set. */ }; /* end of the function checkIdentification. */ getSingularIdentifications(r, s, vectorOfIdentifications) = { /* Check which of the known identification matrices [a,b; c,d] between non-singular vertices */ /* transport the singular point (z,0) into the singular point (Zeta, 0)*/ local(z, Zeta, g, gz, a,b,c,d, vertexTransport, Numerator, Denominator); vertexTransport = listcreate(44); z = component( eval(totalPointSet[r]),1); Zeta = component( eval(totalPointSet[s]),1); for( k = 1, length( vectorOfIdentifications), g = vectorOfIdentifications[k]; a = g[1,1] ; b = g[1,2] ; c = g[2,1] ; d = g[2,2]; /* Check that gz := (az -b)/(-cz+d) = Zeta. */ /* First make sure that the denominator of gz is nonzero: */ Denominator = -nfbasistoalg(K,nfeltmul(K, c,z)) +nfbasistoalg(K, d); if( Denominator != 0, Numerator = nfbasistoalg(K, nfeltmul(K, a,z)) -nfbasistoalg(K, b); gz = nfeltdiv(K, Numerator, Denominator); if( gz == Zeta, /* Record the matrix g into the transportator set. */ listput(vertexTransport, g); ); ); ); /* return the list */ vertexTransport }; /* end of the function getSingularIdentifications */ IsLoop( PointsOnLine:list) = { /* We check if the segment connecting the first two points on the given line */ /* becomes a loop in the quotient space: */ /* This is the case if the two points are identified by some element of the Bianchi group. */ local( pointOne, pointTwo, stabilizerOne, stabilizerTwo, areIdentified); if( length( PointsOnLine) > 1, pointOne = PointsOnLine[1]; pointTwo = PointsOnLine[2]; if( mIs3mod4, stabilizerOne = computeVertexStabilizer3mod4( component( pointOne,1), component( pointOne,2)); stabilizerTwo = computeVertexStabilizer3mod4( component( pointTwo,1), component( pointTwo,2)); ,/* else not congruent 3 mod 4 */ stabilizerOne = computeVertexStabilizer( component( pointOne,1), component( pointOne,2)); stabilizerTwo = computeVertexStabilizer( component( pointTwo,1), component( pointTwo,2)); ); if( length( stabilizerOne) == length( stabilizerTwo), areIdentified = length( getIdentificationMatrices( pointOne, pointTwo, length( stabilizerOne))); ,/* else they can not be identified */ areIdentified = 0; ); , /* else there are no points to identify */ areIdentified = 0; ); /* Return */ areIdentified }; /* end of function IsLoop */ IsHorizontal( Line) = { /* It suffices to check if the imaginary part of the direction vector vanishes. */ if( component( component( Line, 2), 2) == 0, /* return answer: true */ 1 , /* else return answer: false */ 0 ) }; /* end of diagnosis function IsHorizontal */ getSetOfPointsOfSphere(j) = { /* Take only the first component and omit the height. */ local( setOfPoints:list); setOfPoints = listcreate( length( pointsOfSphere[j])); for( r = 1, length( pointsOfSphere[j]), listput( setOfPoints, component( pointsOfSphere[j][r], 1)); ); /* return */ Set( setOfPoints) }; /* end of function getSetOfPointsOfSphere */ identifyIdeals( transportCardinality) = /* Identify singular points which represent the same ideal class. */ { local( z, k, denominatorOf_z, idealAssociatedToCusp); idealAssociatedToCusp = vector( length(totalPointSet)); /* Record into this vector a matrix specifying the associated ideal for every singular point, */ /* and a zero for the points of totalPointSet which are inside Hyperbolic Three-Space. */ for( j = 1, length(totalPointSet), /* Check if point number j lies on the Riemann Sphere: */ if( component(eval(totalPointSet[j]),2) == 0, /* Then call "z" the position of point number j. */ z = component(eval(totalPointSet[j]),1); /* Get the entry of the list singularDenominator, whose first component is z. */ k = 1; while( singularDenominator[k][1] != z, if( k == length(singularDenominator), print("***Error in function identifyIdeals: Pseudo-singular point in fundamental domain."); ); k++; ); denominatorOf_z = singularDenominator[k][2]; /* Record the ideal generated by denominator and numerator of z: */ idealAssociatedToCusp[j] = idealadd(K, denominatorOf_z, z*denominatorOf_z ); ); ); /* Add an identification where the associated ideal is the same: */ for( j = 1, length(totalPointSet), /* and not zero (points inside Hyperbolic Three-Space): */ if ( idealAssociatedToCusp[j] != 0, for( k = j, length(totalPointSet), if ( idealAssociatedToCusp[j] == idealAssociatedToCusp[k] , transportCardinality[j,k] ++; transportCardinality[k,j] ++; ); ); ); ); /* Return the added identifications */ transportCardinality }; /* end of function identifyIdeals */ checkWhichPointsAreInStripe( oldBarycenter, Barycenter_k) = { local( inStripe, RealPart_k, oldRealPart); /* Check which of the two given points is in the preferred stripe of width 1/2. */ if( mIs3mod4 == 0, /* m not congruent to 3 mod 4 */ RealPart_k = component( Barycenter_k, 1); oldRealPart = component( oldBarycenter, 1); if( RealPart_k > 1/2 && oldRealPart <= 1/2, inStripe = "old" ); if( RealPart_k > 1/2 && oldRealPart > 1/2, inStripe = "none" ); if( RealPart_k <= 1/2 && oldRealPart <= 1/2, inStripe = "both" ); if( RealPart_k <= 1/2 && oldRealPart > 1/2, inStripe = "k" ); , /* else, for m congruent to 3 mod 4 */ RealPart_k = component( Barycenter_k, 1) -component( Barycenter_k, 2)/2; oldRealPart = component( oldBarycenter, 1) -component( oldBarycenter, 2)/2; if( RealPart_k >= 0 && oldRealPart < 0 , inStripe = "k" ); if( RealPart_k >= 0 && oldRealPart >= 0 , inStripe = "both" ); if( RealPart_k < 0 && oldRealPart < 0 , inStripe = "none" ); if( RealPart_k < 0 && oldRealPart >= 0 , inStripe = "old" ); ); /* Return */ inStripe }; /* end of function checkWhichPointsAreInStripe */ intersectionOfIdentifications(j,k,r) = { /* Return the following set. */ setintersect( vertexIdentifications[ cornersOf2cell[j][r], cornersOf2cell[k][r]], vertexIdentifications[ cornersOf2cell[j][r+1], cornersOf2cell[k][r+1]] ) }; /* end of function intersectionOfIdentifications */ checkIdentificationIntersection(j,k) = { /* Check if the two cells number j and k are identified by isometries from Gamma. */ local( zeroCellTransportator, prospectiveCellTransportator, r, SHIFT, emptySet, areIdentified); emptySet = Set(); /* Intersect all the sets of isometries transporting a vertex of two-cell number j to a vertex of two-cell number k.*/ /* The entries in the vectors cornersOf2cell have been sorted by orbit numbers. */ /* As a 2-cell might have two corners which lie on the same orbit, there can be pairs */ /* appearing in this sorted vector. We then have to check both the identification suggested */ /* by the sorting, as well as a crossed identification possibility. */ /* To do the latter, we will swap to entries in the cornersOf2cell ordering. */ r = 1; prospectiveCellTransportator = intersectionOfIdentifications(j,k,r); if( prospectiveCellTransportator == emptySet, /* This means that there is no identification, except */ /* if we can get one by swapping the corner ordering. */ swapCorners(k,r,r+1); prospectiveCellTransportator = intersectionOfIdentifications(j,k,r); if( prospectiveCellTransportator == emptySet, /* Last try by swapping back, and then swapping the following pair. */ swapCorners(k,r,r+1); swapCorners(k,r+1,r+2); prospectiveCellTransportator = intersectionOfIdentifications(j,k,r); ); ); zeroCellTransportator = prospectiveCellTransportator; while( r < length( cornersOf2cell[j]) && zeroCellTransportator != emptySet, prospectiveCellTransportator = setintersect( zeroCellTransportator, vertexIdentifications[ cornersOf2cell[j][r+1], cornersOf2cell[k][r+1]]); if( prospectiveCellTransportator == emptySet, /* This means that there is no common identification, except */ /* if we can find another corner on the same orbit, for which this is not empty. */ swapCorners(k,r,r+1); prospectiveCellTransportator = setintersect( zeroCellTransportator, vertexIdentifications[ cornersOf2cell[j][r+1], cornersOf2cell[k][r+1]]); if( prospectiveCellTransportator == emptySet && r+2 <= length( cornersOf2cell[j]), /* Last try by swapping back, and then swapping the following pair. */ swapCorners(k,r,r+1); swapCorners(k,r+1,r+2); prospectiveCellTransportator = setintersect( zeroCellTransportator, vertexIdentifications[ cornersOf2cell[j][r+1], cornersOf2cell[k][r+1]]); ); ); zeroCellTransportator = prospectiveCellTransportator; r++; ); if ( zeroCellTransportator != emptySet, /* the 2-cells number j and k */ areIdentified = 1; ,/* else */ areIdentified = 0; ); /* Return */ areIdentified }; /* end of function checkIdentificationIntersection */ checkBarycenterIdentification(j,k) = { /* Check if there exists an identification between the */ /* barycenters (in the sense described below) of the 2-cells number j and k. */ local( barycenter_j, barycenter_k, areIdentified); /* Construct the barycenters of the projections of the vertices on the 2-cell */ /* onto the complex plane. */ barycenter_j = getBarycenter(j); barycenter_k = getBarycenter(k); /* Find the height of the hemisphere twoCellSupport[j] over barycenter_j. */ heightsquare_j = (1 -norm( nfbasistoalg(K, nfeltmul(K, Mu[twoCellSupport[j]], barycenter_j) -Lambda[twoCellSupport[j]]))) /norm(nfbasistoalg(K,Mu[twoCellSupport[j]])); /* Find the height of the hemisphere twoCellSupport[k] over barycenter_k. */ heightsquare_k = (1 -norm( nfbasistoalg(K, nfeltmul(K, Mu[twoCellSupport[k]], barycenter_k) -Lambda[twoCellSupport[k]]))) /norm(nfbasistoalg(K,Mu[twoCellSupport[k]])); if( mIs3mod4, areIdentified = length( getIdentificationMatrices3mod4( [barycenter_j,heightsquare_j], [barycenter_k,heightsquare_k], 2) ); ,/* else */ areIdentified = length( getIdentificationMatrices( [barycenter_j,heightsquare_j], [barycenter_k,heightsquare_k], 2) ); ); if( areIdentified > 0, areIdentified = 1); /* Return that the cells number j and k */ areIdentified }; /* end of function checkBarycenterIdentifications */ getPreBarycenter(j) = { local(preBarycenter); /* Construct the barycenter of the projections of the vertices on the hemisphere number j onto the complex plane. */ preBarycenter = 0; for ( r = 1, length( pointsOfSphere[j]), preBarycenter = preBarycenter + component( pointsOfSphere[j][r], 1); ); preBarycenter = preBarycenter / length( pointsOfSphere[j]); /* return */ preBarycenter }; /* end of function getPreBarycenter */