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) );
if(maxnormb < normb,maxnormb=normb);
print(maxnormb)
}

trap()
counter = 1
cancel = 0
/* maxnormb is the maximum value of normb from templatemaxnormb */
maxnormb = 4.91465925345340672832447506540745143307521369721 E-27
/* nummflds is number of manifolds in cusped census */
/* nummflds = 4929 */ 
{ for(counter=1,9858,
/* print r c manifold identification ( 5 or 6 or 7 plus number)*/
	ID = counter;
/* print "pr sol" */
/*      sol = extern("sed -n "counter"p perlout1");*/
/*	if(sol == "solutiontype:geometric", ,cancel=1 && print(ID," not geometric"));*/
/*	counter = counter+1; */
/* print "pr sh" */
	va = Vec(extern("sed -n "counter"p perlout1"));
      counter=counter+1;
/* once get shape vector, find n=length of vector */
	vn = matsize(va)[2];
	soltype=2;
      for(imm=1,vn,if(imag(va[imm])<0,soltype=3));
      for(imm=1,vn,if(imag(va[imm])<10^(-28),soltype=4));
/* print "pr fill" */
	FG = extern("sed -n "counter"p perlout1");
/* find total number of equations: n+k */
	numalleq = matsize(FG)[1];
/* find number of cusps */
	vk = numalleq - vn;
/* initialize the cusp equations matrix */
	vF = matrix(vk,2*(vn) +1,ii,jj,FG[ii,jj]);
/* initialize the consistency equation matrix */ 
	vG = matrix(vn,2*(vn) +1,ii,jj,FG[vk+ii,jj]);	
/* if manifold not geometric solution, do not call function, otherwise call funtion */
/*	if(cancel,cancel=0,template(vn,vk,vF,vG));*/
	if( (counter/2)<= 415,
		tet="5 ";
		mfld=(counter/2)-1,
		newcounter=counter - 830;
		if( (newcounter/2)<=962,
			tet="6 ";
			mfld=(newcounter/2)-1,
			tet="7 ";
			mfld=(newcounter/2) - 963
		)
	);
	if(soltype<3,template(vn,vk,va,vF,vG,tet,mfld),
		print("manifold ",tet,mfld," has param with neg or <10^(-28) imag component as approx sol")
	);
    )
}
print("The maximum value for normb for the entire census is ",maxnormb)