v(rr) = vector( (2*n)+1,l,G[rr,l] ) t(ss) = vector( 2*n,l,K[ss,l] ) f(ii) = sum( j = 1, n, F[ii,j]*log(a[j]) ) + sum( j = 1, n, F[ii,n+j]*log(1-a[j]) ) + F[ii,(2*n)+1]*Pi*I g(ii,jj) = ( F[ii,jj]/a[jj] ) - ( F[ii,n+jj]/(1-a[jj]) ) delta(ii,jj) = if( abs(F[ii,jj])+abs(F[ii,jj+n]), ( Eps/(2*n) )*( 1/( (abs(F[ii,jj])/(abs(a[jj]))^2) + \ (abs(F[ii,n+jj])/(abs(1-a[jj]))^2) ) ),((4*normb)/Eta)+1 ) vecdelta(jj) = vector( n, i, delta(i,jj) ) mindelta(jj) = vecmin( vecdelta(jj) ) vecmindelta(jj) = [mindelta(jj),imag(a[jj]),(norml2(a[jj]))/2,(norml2(1-a[jj]))/2] mmindelta(jj) = vecmin( vecmindelta(jj) ) c(i,j) = (abs(F[i,j])/(abs(a[j]) – 2*normhh)^2) + (abs(F[i,j+n])/(abs(1-a[j]) – 2*normhh)^2) template(n,k,a,FF,GG,ttet,mmfld,compare=5) = {local(F,G,H,K,r); F=FF; G=GG; /* read the file FILENAME into SNAP */ /* see that there are h unsurgered cusps */ /* (,) (,) (,) (,) */ /* print shapes - the triangulation has n tetrahedra */ /* print filling equations. Count the number of equations, subtract n from it, and set k equal \ to this difference. First the k-h cusp surgery equations will appear, followed by the h \ meridienal completeness equations, and finally all of the n consistency equations. If no \ surgery has been performed on the manifold, the first k equations are all meridienal \ completeness equations.*/ /* initialize the matrix F, using the first k equations (completeness and surgery equations) */ /* if F is a vector, change it to a one row matrix by \ F = Mat( [X1,1, ..., X1,2n+1] ) */ /* define matrix H, by eliminating the last column of F representing the Pi*I coefficient */ H = matrix (k,2*n,i,j,F[i,j]); /* define the matrix G, using the n consistency equations and\ define matrix K, by eliminating the last column of G representing the Pi*I coefficient */ K = matrix (n,2*n,i,j,G[i,j]); /* redefine F and H by adding rows to them from G and K respectively\ until the rank of F and H are both n */ r=1; /* v(r)= vector( (2*n)+1,l,G[r,l] );*/ /* t(r)=vector( 2*n,l,K[r,l] );*/ while( n – matrank(H) && (n+1-r), if( (matrank(concat(F,v(r))) – matrank(F)) && (matrank(concat(H,t(r))) – matrank(H)), F = concat(F,v(r)); H = concat(H,t(r)), r=r++ ) ); eval(F); eval(H); /* set up the filling equations as log functions evaluated at a*/ /* f(i) = sum( j = 1, n, F[i,j]*log(a[j]) ) + sum( j = 1, n, F[i,n+j]*log(1-a[j]) ) + F[i,(2*n)+1]*Pi*I; */ /* define the vector b in C^n */ b = vector( n, i, f(i) ); /* identify the norm of b */ normb = sqrt( norml2(b) ); /* identify A, the derivative matrix for f at a */ /* g(i,j) = ( F[i,j]/a[j] ) - ( F[i,n+j]/(1-a[j]) ); */ A = matrix( n, n, i, j, g(i,j) ); /* check that determinant of A is not zero */ matdet(A); /* find eigen values for M=(conjugate of A) * (transpose of A) */ M = conj(A)*mattranspose(A); wapprox = polroots( charpoly(M,x) ); w = real( wapprox ); /* define Eta = (square root of smallest eigen value of M)/2 */ Eta = sqrt ( vecmin( w ) )/2; /* define Eps = ((square root of smallest eigen value of M) - Eta)/n */ Eps = ( sqrt( vecmin( w ) ) - Eta )/n; /* Identify DELTA */ /* If not both F[i,j ] and F[i,j+n] are equal to zero, select delta(i,j) according to instructions.\ If both equal zero, then we get zero as the denominator; however any arbitrarily large delta will \ work for this i and j, so select delta(i,j)=(4*normb)/Eta + 1 */ /* delta(i,j) = if( abs(F[i,j])+abs(F[i,j+n]), ( Eps/(2*n) )*( 1/( (abs(F[i,j])/(abs(a[j]))^2) + \ (abs(F[i,n+j])/(abs(1-a[j]))^2) ) ),((4*normb)/Eta)+1 ); */ /* vecdelta(j) = vector( n, i, delta(i,j) ); */ /* mindelta(j) = vecmin( vecdelta(j) ); */ /* vecmindelta(j) = [mindelta(j),imag(a[j]),(norml2(a[j]))/2,(norml2(1-a[j]))/2]; */ /* mmindelta(j) = vecmin( vecmindelta(j) ); */ VECDELTA = vector( n, j, mmindelta(j) ); DELTA = vecmin( VECDELTA ); /* find the value that the norm of b must be less than */ (Eta*DELTA)/4; /* change b into a matrix to do matrix multiplication */ B = matrix(n,1,j,i,b[j]); /* define the vector hh */ hhh = -(A)^(-1)*(B); hh = vector(n, j, hhh[j,1]); /* compare length of hh to DELTA to see if Kantorovich region is included in the Whitney region */ normhh =sqrt(norml2(hh)); 2*(normhh) < DELTA; /* perform the first three tests to see if this method is applicable */ atilde = a + hh; /* test 1 to see if fat solution; if j > n */ for(j = 1, n, if(normhh < imag(atilde[j]), ,error(“failure at atilde[“, j, “]”))); /* test 2 to see if cijj can be defined */ for(j = 1, n, if(normhh < (1/2)*abs(a[j]), ,error(“failure at atilde[“, j, “]”))); /* test 3; other test to see if cijj can be defined */ for(j = 1, n, if(normhh < (1/2)*abs(1 - a[j]), ,error(“failure at atilde[“, j, “]”))); /* identify the Lipschitz ratio, Mlips */ /* c(i,j) = (abs(F[i,j])/(abs(a[j]) – 2*normhh)^2) + (abs(F[i,j+n])/(abs(1-a[j]) – 2*normhh)^2); */ Mlips = sqrt( sum( j = 1, n, sum(i = 1, n, c(i,j)^2) ) ); /* identify normAinv, the norm of A^(-1), using the definition of matrix norm as the supremum of \ A^(-1)v for v on the n-sphere */ normAinv = 1/(2 * Eta); /* do the Kantorovich tests */ /* find the value that the norm of b must be less than or equal to wrt supremum norm */ 1/(2*(normAinv)^2 * Mlips ); /* find the length norm and the value that the norm of b must be less than or equal \ to wrt length norm */ sqrt(norml2(A^(-1))); 1/(2* norml2(A^(-1)) * Mlips ); veccomp=[(Eta*DELTA)/4,1/(2*(normAinv)^2 * Mlips ),1/(2* norml2(A^(-1)) * Mlips )]; maxcomp=vecmax(veccomp); if(maxcomp