/* set precision to 60 (or higher for very large manifolds) from the default of 28 */
\p 60

/* read the file FILENAME into SNAP */
/* see that there are h unsurgered cusps */
/* (,) (,) (,) (,) */

/* print shapes - the triangulation has n tetrahedra */
/* enter the shapes as a vector, so it is regarded as a (1Xn) matrix */
a =
[a1,…,an]
/* find number of tetrahedra */
n = matsize(a)[2]

/* print  filling equations and use this to initialize the matrix FG.  The first k-h equations \
are cusp surgery equations, followed by the h meridianal completeness equations, and finally \
all of the n consistency equations.   */
FG =
[X11,...,X1(2n+1);...;X(n+k)1,...,X(n+k)(2n+1)] 

/* find number of equations (n+k) derived from the "pr fill" command in Snap*/
numalleq = matsize(FG)[1]
/* find total number of cusps, k */
k = numalleq - n

/* initialize F, the cusp equations matrix */
F = matrix(k,2*(n) +1,ii,jj,FG[ii,jj])

/* initialize G, the consistency equation matrix */ 
G = matrix(n,2*(n) +1,ii,jj,FG[k+ii,jj])

/* 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 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+1))
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 D=(conjugate of A) * (transpose of A) */
D = conj(A)*mattranspose(A)
wapprox = polroots( charpoly(D,x) )
w = real( wapprox ) 

/* define Eta = (square root of smallest eigen value of D)/2 */
Eta = sqrt ( vecmin( w ) )/2

/* define  Eps = ((square root of smallest eigen value of D) - 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

/* compare norm of b to above */
normb < ( (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, Lips */
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)
Lips = 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 * Lips )
normb <= 1/(2*(normAinv)^2 * Lips )

/* 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)) * Lips )
(normb  < 1/(2* norml2(A^(-1)) * Lips )) | (normb  = 1/(2* norml2(A^(-1)) * Lips ))