/*
 * The code TQFT is based on the recoupling formulae of the paper:
 * Gregor Masbaum & Pierre Vogel, "$3$-valent graphs and the Kauffman bracket",
 * Pacific J. Math.  164  (1994),  no. 2, 361--381,
 * and  computes the natural action of the
 * mapping class group on the so(su)-Verlinde modules 
 * V_k(S) of level k of a surface S.
 * 
 * The code uses Pari-gp, in particular the very 
 * efficient computations in number
 * fields. Before running TQFT, please install Pari on your system. 
 *
 * The following functions are used in the code:
 * init_so(k,boolean), init_su(k) initialize the global variable 
 * KLEVEL and TYPE: in case su KLEVEL=4*k+8 and TYPE=-k 
 * or in case so TYPE=k and KLEVEL=2*k+4 with k odd,
 * and moreover the global variable POL="KLEVEL"-cyclotomic polynomial. All
 * computations are done in the number field NF=Q[A]/(POL). This corresponds
 * to the specialization of the quantum variable A of the Kauffman bracket 
 * to a primitive KLEVEL'th root of unity. 
 * If boolean=1, the code will work with matrices in sparse format. The
 * default is boolean=0.
 * Further functions are: qi(n,poly) (and qi_alt(n,poly)), qfact(n,poly)
 * mu(i,poly), tomega_su(poly), tomega_so(poly). 
 * Important are the coefficient functions
 * triheco(a,b,c,poly), tetraheco(a,b,c,d,e,f,poly), see the paper 
 * of Masbaum and Vogel.
 * The coefficient functions triheco(a,b,c,poly) and 
 * tetraheco(a,b,c,d,e,f,poly) 
 * give an answer  only  if at the triple points
 * of the graph the three incident "colors" a,b,c satisfy 
 * conditions, that are tested
 * by the function admiss(a,b,c,k).
 * The dimensions grow exponentially with the level k. So, we implemented also
 * a sparse version. A related global variables is  DO_S_matr. The function 
 * mult_sparse(t,s) multiplies the tabular matrix t and sparse matrix s. 
 * Given a rooted
 * planar tree by a vector, say [0,1,1,3,3], meaning that the first vertex
 * is attached to the root(=0), the next two to the first(=1), and last two
 * to the third(=3). The function init_boom([0, ...],inp) initializes
 * the labelling of the corresponding trivalent graph, where the free edge
 * is colored by the color inp. The function init_boom(boom,inp)
 * initializes also the global variable BoomBasis, which is a list of
 * admissible edge colorings. See the pictures in tqft.pdf.
 */

\\ INITIALIZE KLEVEL AND NUMBERFIELD Q[A]/(POL). 
\\ Think of "A" as "Kauffman's A"
\\ and POL as an appropriate cyclotomic polynomial.

global(DO_S_matr,KLEVEL,A,POL,NF,BOOM,TV,TE,TOM_so,TOM_su,TYPE=0,DO_Retrieve);
init(k,do_s=0)=
{
if(k>0, init_so(k , do_s), init_su(-k , do_s))
}
init_so(k , do_s=0) = \\ k odd integer, do_s boolean.
{
       local(str);
       if( divrem(k,2)[2], DO_S_matr=do_s;
			   if( TYPE != k, 
			       KLEVEL=k;POL=polcyclo(2*KLEVEL+4,A);
			       NF=bnfinit(POL);
			         if( KLEVEL > 19 ,DO_Retrieve=0,
				DO_Retrieve=1;
				str=concat(["zcat tt/tetrah_so_",KLEVEL,".gz"]);
				TV=extern(str);	    
				str=concat(["zcat tt/trih_so_",KLEVEL,".gz"]);
				  TE=extern(str));
				 TOM_so=tomega_so(POL);
				 );
			   TYPE=k;     \\ if TYPE > 0, we do "so".
			   print("so, "," KLEVEL= ",KLEVEL,", POL= ",POL),
       print("WARNING: k  not odd")     )
}
init_su(k, do_s = 0 ) =       \\ k integer, k > 1, do_s boolean.
{
      local(str);
      DO_S_matr=do_s;
      if( TYPE != -k,
            KLEVEL=k;POL=polcyclo(4*KLEVEL+8,A);
	    NF=bnfinit(POL);
	    if( KLEVEL > 14 ,DO_Retrieve=0,
	    DO_Retrieve=1;
	    str=concat(["zcat tt/tetrah_su_",KLEVEL,".gz"]);
	    TV=extern(str);
	    str=concat(["zcat tt/trih_su_",KLEVEL,".gz"]);
	    TE=extern(str));
	    TOM_su=tomega_su(POL);
	    );
      TYPE=-k;	             \\ if TYPE < 0, we do "su".
      print("su, "," KLEVEL= ",KLEVEL,", POL= ",POL)
}
Verlinde(g,k)=     \\ Dimension of the so-module for genus g, KLEVEL=k,
                   \\ without punctures.
{
local(b);
if(k<0,return("Only so, but su will come"));
pol=polcyclo(2*k+4,a);q=Mod(a^2,pol);
lift((-(k+2)/4)^(g-1)*sum(j=1,(k+1)/2,((q^j-1/q^j)/2)^(2-2*g)))
}
VerlindePol(g)=   \\ computes the Verlinde_so dimension as polynomial
                  \\ in the number h of even colors: 0,2, ... ,2*(h-1).
		  \\ h=(KLEVEL+1)/2. The output is the Ehrhardt polynomial
		  \\ of the polytope defined by the inequalities for
		  \\ admissibility.
{
local(X,Y);
X=vector(3*g,h,h);Y=vector(3*g,h,Verlinde(g,2*h-1));
polinterpolate(X,Y,h)
} 

\\ USEFULL "QUANTUM" QUANTITIES IN THE NUMBER FIELD Q[A]/(poly)

qi(n,poly)=         \\ quantum integer n, 0 <= n <= k. 
{
	Mod(A^(2*n)-A^(-2*n),poly)/Mod(A^2-A^(-2),poly);
}
qi_alt(n,poly)=
{
       (-1)^n*qi(n+1,poly)
}
qfact(n,poly)=      \\ quantum factoriel.
{
        prod(i=1,n,qi(i,poly))
}
mu(i,poly)=
{
        Mod(-A,poly)^(i^2+2*i)
}
tomega_so(poly)=
{
        sum(i=0,(KLEVEL-1)/2,mu(2*i,poly)*qi(2*i+1,poly)^2)
}
tomega_su(poly)=
{
        sum(i=0,KLEVEL,mu(i,poly)*qi(i+1,poly)^2)
}

\\ 3- AND 6- SYMBOLS AND COEFFICIENTS IN THE NUMBER FIELD Q[A]/(poly)

triheco(a,b,c,poly)=        \\ a,b,c admissible tripel of colors
{
        local(s,ss,aa,bb,cc,k);
        k=KLEVEL; s=a+b+c;
        if( admiss(a,b,c,k) ,
        if( a*b*c > 0,
          aa=(s-2*a)/2;bb=(s-2*b)/2;cc=(s-2*c)/2;ss=aa+bb+cc;
        (-1)^ss*qfact(ss+1,poly)*\
	qfact(aa,poly)/qfact(bb+cc,poly)*\
	qfact(bb,poly)/qfact(aa+cc,poly)*\
	qfact(cc,poly)/qfact(aa+bb,poly),
        (-1)^(s/2)*qi(s/2+1,poly)),     0      )

}
hatwco(a,b,c,poly)=
{
        local(e,s,aa,bb,cc,k);
        k=KLEVEL;s=a+b+c;
        aa=(s-2*a)/2; bb=(s-2*b)/2; cc=(s-2*c)/2; ss=aa+bb+cc;
        e=aa*bb-cc*(ss+2);
       (-1)^cc*Mod(A,poly)^e
}
linking(i,j,poly)=
{
(-1)^(i+j)*Mod((A^(2*(i+1)*(j+1))-A^(-2*(i+1)*(j+1)))/(A^2-1/A^2),poly)
}
admiss(a,b,c,k)=      \\ boolean test function for admissibility.
{
       local(aa,bb,cc,s);
       s=a+b+c; 
       if(s <= 2*k && divrem(s,2)[2] == 0,
       aa=(s-2*a)/2; bb=(s-2*b)/2; cc=(s-2*c)/2;
       if(0 <= aa && 0 <= bb && 0 <= cc , 1 , 0 ) , 0)
}
trih_so()=       \\ precomputes values for triheco, see init_so(k).
{
        local(S,m,a,b,c,z,N,k,h);
	k=KLEVEL;poly=polcyclo(2*k+4,A);
	h=(k+1)/2;
        bas=vector(h^3);
	trih=vector(h^3);
                N=3;
        for(s=1,h^3,
                S=s-1;m=vector(N);
                for(n=1,N,
                        m[n]=divrem(S,h^(N-n))[1];
                        S=S-m[n]*h^(N-n);
                        );
                 m=2*m;a=m[1];b=m[2];c=m[3];
       if(admiss(a,b,c,k) , z++;
                     trih[z]=[s,m,triheco(a,b,c,poly)];);
       );
       trih=vector(z,i,trih[i]);
}
trih_su()=      \\ precomputes values for triheco, see init_su(k)
{
        local(S,m,a,b,c,z,N,k,h,poly);
        k=KLEVEL;poly=polcyclo(4*k+8,A);
        h=k+1;
        trih=vector(h^3);
                N=3;
        for(s=1,h^3,
                S=s-1;m=vector(N);
                for(n=1,N,
                        m[n]=divrem(S,h^(N-n))[1];
                        S=S-m[n]*h^(N-n);
                        );
                 a=m[1];b=m[2];c=m[3];
       if(admiss(a,b,c,k) , z++;
                      trih[z]=[s,m,triheco(a,b,c,poly)];);
       );
       trih=vector(z,i,trih[i]);
}

tetrah_so()=    \\ precomputes values of tetraheco, see init_so().
{     
        local(S,m,a,b,c,d,e,f,z,N,h,poly);
	poly=polcyclo(2*KLEVEL+4,A);
        h=(KLEVEL+1)/2;
	tetrah=vector(h^6);
                N=6;
        for(s=1,h^6,

                S=s-1;
                m=vector(N);
                for(n=1,N,
                        m[n]=divrem(S,h^(N-n))[1];
                        S=S-m[n]*h^(N-n);
                        );
                 m=2*m;
                 a=m[1];b=m[2];e=m[3];d=m[4];c=m[5];f=m[6];
        if(admiss(a,b,e,KLEVEL) && admiss(b,d,f,KLEVEL) && \
                admiss(c,d,e,KLEVEL) && admiss(a,c,f,KLEVEL) ,
        z++;
	tetrah[z]=[s,m,tetraheco(a,b,e,d,c,f,poly)])
	   );
	tetrah=vector(z,i,tetrah[i]);
}
tetrah_su()=    \\ precomputes values of tetraheco, see init_su().
{
        local(S,m,a,b,c,d,e,f,z,N,h,poly);
	poly=polcyclo(4*KLEVEL+8,A);
        h=KLEVEL+1;
        tetrah=vector(h^6);
                N=6;
        for(s=1,h^6,

                S=s-1;
                m=vector(N);
                for(n=1,N,
                        m[n]=divrem(S,h^(N-n))[1];
                        S=S-m[n]*h^(N-n);
                        );
                
                 a=m[1];b=m[2];e=m[3];d=m[4];c=m[5];f=m[6];
        if(admiss(a,b,e,KLEVEL) && admiss(b,d,f,KLEVEL) && \
                admiss(c,d,e,KLEVEL) && admiss(a,c,f,KLEVEL) ,
        z++;
        tetrah[z]=[s,m,tetraheco(a,b,e,d,c,f,poly)])
           );
        tetrah=vector(z,i,tetrah[i]);
}

tetraheco(a,b,e,d,c,f,poly)=     \\ complicated coefficient.
{
       local(abe,bfd,cde,acf,ad,ef,bc,s,MAX,MIN,FACTOR,SOM,k,m);
           k=KLEVEL;
       if( admiss(a,b,e,k) && admiss(b,d,f,k) && \
	       admiss(c,d,e,k) && admiss(a,c,f,k) ,
           s=a+b+c+d+e+f;
           abe=(a+b+e)/2; bfd=(b+f+d)/2;
           cde=(c+d+e)/2; acf=(a+c+f)/2;
           ad=(s-a-d)/2; ef=(s-e-f)/2; bc=(s-b-c)/2;
             MAX=max(max(abe,bfd),max(cde,acf));
             MIN=min(ad,ef);MIN=min(MIN,bc);
	     SOM=sum(m=MAX,MIN,(-1)^m*qfact(m+1,poly)\
  /qfact(ad-m,poly)/qfact(ef-m,poly)/qfact(bc-m,poly)\
 /qfact(m-abe,poly)/qfact(m-bfd,poly)/qfact(m-cde,poly)/qfact(m-acf,poly));
      if( SOM == 0 ,   0 ,    
         FACTOR=1/qfact(a,poly)/qfact(b,poly)/qfact(c,poly)\
	     /qfact(d,poly)/qfact(e,poly)/qfact(f,poly)\
   *qfact(ad-abe,poly)*qfact(ef-abe,poly)*qfact(bc-abe,poly)\
   *qfact(ad-cde,poly)*qfact(ef-cde,poly)*qfact(bc-cde,poly)\
   *qfact(ad-acf,poly)*qfact(ef-acf,poly)*qfact(bc-acf,poly)\
   *qfact(ad-bfd,poly)*qfact(ef-bfd,poly)*qfact(bc-bfd,poly);
   FACTOR*SOM), print(["WARNING: not admissable:",[a,b,e,d,c,f]])   )
}
zes_j_symb(a,b,i,c,d,j,poly)=    \\ 6_j_symbol. zes=6 in Dutch. 
{
      qi_alt(i,poly)/rtriheco(i,a,d,poly)/rtriheco(i,b,c,poly)*\
      rtetraheco(i,b,c,j,d,a,poly)
}
rtetraheco(a,b,e,d,c,f,poly)=  \\If DO_Retrieve==1, Retrieves precomputed 
                                \\tetraheco  values from  TV in database tt.
{
local(z,m,h,j,u,v);
if(DO_Retrieve==0,return(tetraheco(a,b,e,d,c,f,poly)));
if(TYPE < 0,
h=KLEVEL+1;
m=[a,b,e,d,c,f],
h=(KLEVEL+1)/2;
m=[a/2,b/2,e/2,d/2,c/2,f/2]);
z=1+sum(i=1,6,m[i]*h^(6-i));
u=1;v=matsize(TV)[2];j=1;
while(v-u > 1, j=floor((u+v)/2);if(TV[j][1]<z, u=j, v=j));
if(TV[u][1]==z, j=u, j=v);m=TV[j][2];
if(m[1]==a && m[2]==b && m[3]==e && m[4]==d && m[5]==c && m[6]==f,   
         TV[j][3]  ,
print(["WARNING:",a,b,e,c,d,f,"NOT ADMISSIBLE"])) 
}
rtriheco(a,b,c,poly)=   \\ If DO_Retrieves==1, retrieves precomputed triheco 
                        \\ values from TE in database tt.
{
local(z,m,h,j,u,v);
if(DO_Retrieve==0,return(triheco(a,b,c,poly)));
if(TYPE < 0,
h=KLEVEL+1;
m=[a,b,c],
h=(KLEVEL+1)/2;
m=[a/2,b/2,c/2]);
z=1+sum(i=1,3,m[i]*h^(3-i));
u=1;v=matsize(TE)[2];j=1;
while(v-u > 1, j=floor((u+v)/2);if(TE[j][1]<z, u=j, v=j));
if(TE[u][1]==z, j=u, j=v);m=TE[j][2];
if(m[1]==a && m[2]==b && m[3]==c , TE[j][3]  ,
print(["WARNING:",a,b,c,"NOT ADMISSIBLE"]))
}
\\ TREE=BOOM(in Dutch)
/*
 * Given a rooted planar tree, we attach to it a surface. For a planar tree
 * we use a code like [0,1,2,2,2]. The function BV(boom) computes a labeling
 * of the vertices, edges and regions of a special doubling of the planar 
 * tree. The doubling is a trivalent graph as in the paper of Masbaum and Vogel
 * The function init_boom(boom,inp) computes a basis of the space of
 * linear combinations with integral coefficients of admissible colorings 
 * of this trivalent graph with colors c in the range
 * 0<= c <=k, and moreover  c even in the so case.
 */

BV(boom)=        \\ boom as a vector [0,1, ... ,]
{
        local(prev,BVprev,res,res1,res2,bL,attach,i_in,i_uit,a,b);
        bL=matsize(boom)[2];
	res1=[1,1,2,1,1,0];
        res2=[3,1,4,1,1,0, 1,2,5,1,2,3, 2,3,5,2,1,2];

        if( bL == 1 , res=matrix(1,6,i,j,res1[j]));
        if( bL == 2  , res=matrix(3,6,i,j,res2[6*(i-1)+j]));
        if( bL > 2,
        prev=vector(bL-1,u,boom[u]);

        BVprev=BV(prev);
        attach=boom[bL];
        if( attach == 1 , i_uit=2*bL-2);
        if( attach > 1 ,    i_in=2;
             while(BVprev[i_in,5]<attach,i_in++);
             i_uit=BVprev[i_in,6];);
             res=matrix(2*bL-1,6);
        for(i=1,i_uit-1,for(j=1,6,res[i,j]=BVprev[i,j]);
                        res[i,3]=res[i,3]+2;
                       if(res[i,6]>i_uit-1,res[i,6]=res[i,6]+2););

        a=vector(6); a=[i_uit-1,i_uit,3*bL-1,attach,bL,i_uit+1];
        for(j=1,6,res[i_uit,j]=a[j]);
        b=vector(6); b=[i_uit,i_uit+1,3*bL-1,bL,attach,i_uit];
        for(j=1,6,res[i_uit+1,j]=b[j]);
        for(i=i_uit+2,2*bL-1,for(j=1,6,res[i,j]=BVprev[i-2,j]);
                             for(j=1,3,res[i,j]=res[i,j]+2);
               if(res[i,6]>i_uit,res[i,6]=res[i,6]+2);
                       );
        res[1,1]=2*bL-1; res[1,3]=2*bL;
        );
        res
}
init_boom(boom,inp)=    \\ boom is a tree code and inp is a color.
{
        local(bL,k,N,r,z,adm,m,BM,BB);
        global(BoomLengte,BoomMatrix,BoomBasis);
	if( TYPE == 0, return("Do first: init_so(level) or init_su(level"));
        k=KLEVEL;if(TYPE>0,kk=(k+1)/2,kk=k+1);
        BoomLengte=matsize(boom)[2];
        bL=BoomLengte;N=3*bL-1;
        BM=BV(boom);BB=vector(kk^(N-1));
        z=0;m=vector(N);
        for(s=1,kk^(N-1),
                S=s-1;
                for(n=1,N-1,
                        if(n<2*bL,r=n,r=n+1);
                        m[r]=divrem(S,kk^(N-n-1))[1];
                        S=S-m[r]*kk^(N-n-1);
                        );
		if(TYPE>0, m=2*m;m[2*bL]=inp , m[2*bL]=inp);
                adm=1;n=1;
                while(adm && n<2*bL,
                        adm=admiss(m[BM[n,1]],m[BM[n,2]],m[BM[n,3]],k);n++);
                if(adm,z++;BB[z]=[m,z]);
                );
        BOOM=1;
        BoomMatrix=BM;
        BoomBasis=vector(z,i,BB[i]);
        z
}
\\ ACTIONS OF DEHN TWISTS
/*
 * The surface to the linear tree [0,1,2, ... ,g-1] is a surface
 * of genus g with one boundary component. The curves b_1,b_2, ... ,b_g
 * surround in the standard picture of the surface the g holes.
 * The function twist_so(su)_boom_B(r) computes the matrix of the
 * positive Dehn twist along b_r on the space of colorings with input color
 * equal to inp as given by init_boom(boom, inp). The function 
 * twist_so(u)_boom_A(n) computes the matrix
 * of the Dehn twists along curves a_n, 1<= n <= 2*g-1, that bound
 * in the handlebody of the standard picture embedded disks in the 
 * handlebody. 
 * tw[A,B,C]() are short cuts for twist_boom_so_[A,B,C]().
 * For the case boom=[0,1,2, ... ,g-1] the Dehn twists along b_r,a_n 
 * generate by a result of Lickorish the mapping class group of the surface of 
 * genus g with one boundary component. Remember that the input color
 * has been fixed in the initialization init_boom(u)(boom,inp) 
 * For more general planar trees, one gets a set of a_n's
 * and b_r's that compute in a handy way the twist of slalom knots.
 * See the paper N. A'Campo, "Slalom knots", Publ. Math. IHES 1998.
 */ 

twist_so_boom_A(n)=      \\ Twist about the curve that links 
                         \\ with the nth edge, 1<= n <=3*bL-1.
{
        local(BB,dim);
        BB=BoomBasis; dim=matsize(BB)[2];
	if(DO_S_matr == 1,
	   vector(dim+2,u,
	   if(u == 1, [dim,dim,"M"],
		   if(u<dim+2,[u-1,u-1,1/mu(BB[u-1][1][n],POL)],[0,0,0])) ),
           matrix(dim,dim,i,j,if( i == j ,1/mu(BB[i][1][n],POL) , 0))) ;
}
twA(n)=
{
local(res);
if(TYPE > 0 ,res=twist_so_boom_A(n));
if(TYPE < 0 ,res=twist_su_boom_A(n));
res
}
twist_su_boom_A(n)=      \\ Twist about the curve that links 
                         \\ with the nth edge, 1<= n <=3*bL-1.
{
        local(BB,dim);
        BB=BoomBasis; dim=matsize(BB)[2];
        if(DO_S_matr == 1,
           vector(dim+2,u,
           if(u == 1, [dim,dim,"M"],
                   if(u<dim+2,[u-1,u-1,1/mu(BB[u-1][1][n],POL)],[0,0,0])) ),
           matrix(dim,dim,i,j,if( i == j ,1/mu(BB[i][1][n],POL) , 0))) ;
}
twB(r)=
{
if(DO_S_matr,  twist_boom_B_s(r) , twist_boom_B_t(r));
}
twist_boom_B_t(r)=   \\ Twist about the curve that surrounds
                      \\ the rth hole, 1  <= r <= genus = BoomLengte.
		      \\ Assumes DO_S_matr==0. Tabular output.
{
        local(bL,BM,BB,TOM,ring_size,i_in,ii_in,ecke,\
	     inp_ind,m_ind,N,m_ext,n_ext,mi,nj,M,MM,Mi,C);
        bL=BoomLengte; BM=BoomMatrix; BB=BoomBasis; 
	if(r == 1 && bL == 1, N=1;inp_ind=vector(N,u,2);m_ind=vector(N,u,1));
        if(r == 1 && bL>1,
                ii_in=1; ring_size=0;
                while(ii_in > 0 , ring_size++;
                        ii_in=BM[if( ii_in == 2*bL-1 , 1, ii_in+1),6];);
                N=2*ring_size-1; inp_ind=vector(N); m_ind=vector(N); ecke=1;
		        for(i=1,N,if(divrem(i,2)[2],
                    if( ecke == 1 ,inp_ind[i]=2*bL, inp_ind[i]=BM[ecke,1]);
                    m_ind[i]=BM[ecke,2]; ecke=ecke+1,
                    inp_ind[i]=BM[ecke,2]; m_ind[i]=BM[ecke,3];
                    ecke=BM[ecke,6];)
                            );
          );
        if(  r > 1 ,
                i_in=1;
        while( BM[i_in,3] != 2*bL-1+r ,i_in++);
        ecke=i_in; ii_in=BM[ecke+1,6];
        if(ii_in == i_in, ring_size=1,
		          ring_size=1;
        while(ii_in != i_in,ii_in=BM[ii_in+1,6];ring_size++););
        N=2*ring_size; inp_ind=vector(N); m_ind=vector(N); ecke=i_in;
        for( i=1 , N , if(divrem(i,2)[2],
	      inp_ind[i]=BM[ecke,1]; m_ind[i]=BM[ecke,2]; ecke=ecke+1, 
	      inp_ind[i]=BM[ecke,2]; m_ind[i]=BM[ecke,3]; ecke=BM[ecke,6];));
	  );
        dim=matsize(BB)[2];
		   M=matrix(dim,dim);
		   if(TYPE>0, TOM=TOM_so, TOM=TOM_su);
        for(i=1,dim, 
	         for(j=1,dim,
                mi=BB[i][1]; nj=BB[j][1]; m_ext=mi; n_ext=nj;
                for(n=1,N,m_ext[m_ind[n]]=0; n_ext[m_ind[n]]=0;);
                if( sum(g=1,3*bL-1,(m_ext[g]-n_ext[g])^2) , 
		    M[i,j]=0 ,
      C=sum(h=0,KLEVEL,
      if(TYPE< 0 ||  divrem(h,2)[2]==0,   
      if( prod(n=1,N, admiss(h,mi[m_ind[n]],nj[m_ind[n]],KLEVEL)) ,
      mu(h,POL)*qi_alt(h,POL)*\
      prod(n=1,N, qi_alt(nj[m_ind[n]],POL)*\
      rtetraheco(mi[inp_ind[n]],mi[m_ind[if(n==1,N,n-1)]],\
         mi[m_ind[n]],h,nj[m_ind[n]],\
         nj[m_ind[if(n==1,N,n-1)]],POL)\
   /rtriheco(h,nj[m_ind[n]],mi[m_ind[n]],POL)\
   /rtriheco(mi[inp_ind[n]],nj[m_ind[if(n==1,N,n-1)]],nj[m_ind[n]],POL))
            ,0),0));

	M[i,j]=C/TOM);
	      );
	      );
M
}
twist_boom_B_s(r)=   \\ Twist about the curve that surrounds
                      \\ the rth hole, 1  <= r <= genus = BoomLengte.
		      \\Assumes DO_S_matr==1. Sparse output.
{
        local(bL,BM,BB,TOM,ring_size,i_in,ii_in,ecke,\
             inp_ind,m_ind,N,m_ext,n_ext,mi,nj,M,MM,Mi,C);
        bL=BoomLengte; BM=BoomMatrix; BB=BoomBasis;
        if(r == 1 && bL == 1, N=1;inp_ind=vector(N,u,2);m_ind=vector(N,u,1));
        if(r == 1 && bL>1,
                ii_in=1; ring_size=0;
                while(ii_in > 0 , ring_size++;
                        ii_in=BM[if( ii_in == 2*bL-1 , 1, ii_in+1),6];);
                N=2*ring_size-1; inp_ind=vector(N); m_ind=vector(N); ecke=1;
                        for(i=1,N,if(divrem(i,2)[2],
                    if( ecke == 1 ,inp_ind[i]=2*bL, inp_ind[i]=BM[ecke,1]);
                    m_ind[i]=BM[ecke,2]; ecke=ecke+1,
                    inp_ind[i]=BM[ecke,2]; m_ind[i]=BM[ecke,3];
                    ecke=BM[ecke,6];)
                            );
          );
        if(  r > 1 ,
                i_in=1;
        while( BM[i_in,3] != 2*bL-1+r ,i_in++);
        ecke=i_in; ii_in=BM[ecke+1,6];
        if(ii_in == i_in, ring_size=1,
                          ring_size=1;
        while(ii_in != i_in,ii_in=BM[ii_in+1,6];ring_size++););
        N=2*ring_size; inp_ind=vector(N); m_ind=vector(N); ecke=i_in;
        for( i=1 , N , if(divrem(i,2)[2],
              inp_ind[i]=BM[ecke,1]; m_ind[i]=BM[ecke,2]; ecke=ecke+1,
              inp_ind[i]=BM[ecke,2]; m_ind[i]=BM[ecke,3]; ecke=BM[ecke,6];));
          );
        dim=matsize(BB)[2];
        if(TYPE>0, TOM=TOM_so, TOM=TOM_su);
                   MM=vector(dim+2);
                   MM[1]=vector(1,i,[dim,dim,"M"]);
                   MM[dim+2]=vector(1,i,[0,0,0]);

        Mi=[];
        for(i=1,dim,
                 for(j=1,dim,
                mi=BB[i][1]; nj=BB[j][1]; m_ext=mi; n_ext=nj;
                for(n=1,N,m_ext[m_ind[n]]=0; n_ext[m_ind[n]]=0;);
                if( sum(g=1,3*bL-1,(m_ext[g]-n_ext[g])^2) == 0,
                 
      C=sum(h=0,KLEVEL, if(TYPE<0 || divrem(h,2)[2]==0,
      if( prod(n=1,N, admiss(h,mi[m_ind[n]],nj[m_ind[n]],KLEVEL)) ,
      mu(h,POL)*qi_alt(h,POL)*\
      prod(n=1,N, qi_alt(nj[m_ind[n]],POL)*\
      rtetraheco(mi[inp_ind[n]],mi[m_ind[if(n==1,N,n-1)]],\
         mi[m_ind[n]],h,nj[m_ind[n]],\
         nj[m_ind[if(n==1,N,n-1)]],POL)\
   /rtriheco(h,nj[m_ind[n]],mi[m_ind[n]],POL)\
   /rtriheco(mi[inp_ind[n]],nj[m_ind[if(n==1,N,n-1)]],nj[m_ind[n]],POL)),0),0));

        if( C == 0 ,     ,  Mi=concat(Mi,[[i,j,C/TOM]]) );
         ););
         MM[i+1]=Mi;
         Mi=[];);
         M=concat(MM);
}
twist_boom_C(r)=   \\ Twist about the curve that cuts off
                      \\ the rth hole, 1  < r <= genus = BoomLengte.
{
        local(bL,BM,BB,N,m_ext,n_ext,mi,nj,M,C,p,q,rr,a,b,c,d);
        bL=BoomLengte; BM=BoomMatrix; BB=BoomBasis;
	if( r == 1, M=twist_so_boom_A(2*bL));
        if( r > 1,
	rr=2*bL+r-1;
        p=1;
	while(BM[p,3] != rr,p++);
	q=BM[p,6];
        dim=matsize(BB)[2];
        if(DO_S_matr == 1 ,
	           MM=vector(dim+2);
                   MM[1]=vector(1,i,[dim,dim,"M"]);
		   MM[dim+2]=vector(1,i,[0,0,0]),
                   M=matrix(dim,dim)  );
        Mi=[];
        for(i=1,dim,
            for(j=1,dim,
                mi=BB[i][1]; nj=BB[j][1]; m_ext=mi; n_ext=nj;
                m_ext[rr]=0; n_ext[rr]=0;
                if( sum(g=1,3*bL-1,(m_ext[g]-n_ext[g])^2) ,
                     if(DO_S_matr,  ,M[i,j]=0) ,
      a=mi[BM[p,1]];b=mi[BM[p,2]];c=mi[BM[q,1]];d=mi[BM[q,2]];
      
      C=sum(h=0,KLEVEL,if(TYPE<0 || divrem(h,2)[2]==0,
      if( admiss(h,a,d,KLEVEL)*admiss(h,b,c,KLEVEL),
      1/mu(h,POL)*\
      zes_j_symb(a,b,h,c,d,mi[rr],POL)*\
      zes_j_symb(b,c,nj[rr],d,a,h,POL))  ,  0) ););
      if(DO_S_matr, Mi=concat(Mi,[[i,j,C]]),M[i,j]=C););
      
      if(DO_S_matr, MM[i+1]=Mi;Mi=[]);
            );    
      if( DO_S_matr == 1, M=concat(MM););
      );
      M
}      
twC(r)=
{
twist_boom_C(r)
}

\\ EXTRA
/*
 * The function slalom(boom,k,inp) computes the action of the monodromy 
 * of the slalom knot
 * in TQFT. The function eval_sl(m,k,inp) computes the action 
 * of the integram 2x2 matrix
 * m with determinant=1 on the TQFT-modules of the once-punctured 
 * torus in level k and
 * with input color inp. The function normk(t) computes for 
 * t in Q(A)/POL the maximal
 * possible absolut value of s(t), s being an embedding of Q(A)/POL in C.
 */
slalom(boom,k,inp)=  \\ monodromy "so" or "su" of slalom knot in level
                     \\ acording to sign of k.
		     \\ The result is the TQFT action of the 
		     \\ monodromy of the fibered slalom knot.
{
local(res,dim,m,slalomA);
	if( TYPE != k, 
	if(k>0,init_so(k,1),init_su(-k,1)));
	DO_S_matr=1;
	init_boom(boom,inp);
	dim=matsize(BoomBasis)[2];
	res=s2t(twB(1));
	if(BoomLengte>1,
	for(r=2,BoomLengte,res=mult_sparse(res,twB(r))););
	slalomA=vector(dim+2);slalomA[1]=[dim,dim,"M"];slalomA[dim+2]=[0,0,0];
        for(u=2,dim+1,
        m=mu(BoomBasis[u-1][1][1],POL)*\
        prod(ii=2,BoomLengte,mu(BoomBasis[u-1][1][2*BoomLengte-1+ii],POL) );
        slalomA[u]=[u-1,u-1,1/m]);
	mult_sparse(res,slalomA);
}
coxeter(boom,k,inp)=  \\ Coxeter element c. It satisfies
                      \\ twB(r)=c*twA(2*bL-1+r)/c, r>1, and
                      \\ twB(1)=c*twA(1)/c.
		      \\ The result is up to conjugacy the
		      \\ TQFT action of the monodromy.
{
       local(res,bL); 
       if(k>0,init_so(k,1),init_su(-k,1)); DO_S_matr=1;
       init_boom(boom,inp);
       bL=BoomLengte;
       res=eval_s([twA(1),twB(1)]);
      for(r=2,bL,res=mult_sparse(res,twA(2*bL-1+r));
                 res=mult_sparse(res,twB(r));
      );
      res
}
saturb(boom,k,inp)=    \\ Computes a lattice, i.e. a Z_K module of maximal
                       \\ in the so-module K^dim, dim=matsize(BoomBasis)[2]
                       \\ K=Q(A)/(POL), that is invariant under the TQFT-action
                       \\ of the mapping class group of the punctured surface.
                       \\ The first component of the output 
                       \\ P=saturb(boom,k,inp)[1] is dim x dim matrix who's
                       \\ collums generate the lattice. It follows for instance,
                       \\ that 1/P*coxeter(boom,k,inp)*P has its entries in Z_K.
{
local(dim,idl,gen,stop,iter,test,Ma,Mb,Mg,diag,index,indexZ);
init(k);init_boom(boom,inp);
dim=matsize(BoomBasis)[2];  
idl=vector(2*dim,i,idealhnf(NF,1));
Ma=matid(dim);vol=1;
gen=vector(4*BoomLengte-1);                                                     
for(i=1,4*BoomLengte-1,
      if(i<3*BoomLengte,gen[i]=twA(i),gen[i]=twB(i-(3*BoomLengte-1)) ));
stop=0;iter=0;indexZ=1;
while(stop<1,iter++;
for(g=1,4*BoomLengte-1,Mg=matalgtobasis(NF,concat([gen[4*BoomLengte-g]*Ma,Ma]));
Mb=nfhnf(NF,[Mg,idl]);                                                          
diag=matbasistoalg(NF,matrix(dim,dim,i,j,
if(i==j, bnfisprincipal(NF,Mb[2][i])[2] , 0)  ) );                              
Ma=matbasistoalg(NF,Mb[1])*diag;
);
test=lift(vol/matdet(Ma));
index=1;
for(j=1,2*k+4,if(gcd(j,2*k+4)==1,index=index*subst(test,A,A^j)));
index=lift(Mod(index,POL));
print("iter=",iter,", ","index = ",index,", ","test = ", test);
if(index==1, stop++, vol=matdet(Ma));
indexZ=indexZ*index;
);
[Ma,Mb,iter,indexZ]
}
eval_s(w)=   \\ computes the product w[1]*w[2]* ... w[length(w)]
                \\ assuming that the factors w[i] are sparse matrices.
{
	local(dim,res,i);
	dim=matsize(BoomBasis)[2];
	res=matid(dim);i=1;
	while(i<length(w)+1,res=mult_sparse(res,w[i]);i++);
	res
}
eval_t(w)=   \\ computes the product w[1]*w[2]* ... w[length(w)]
                \\ assuming that the factors w[i] are tabular matrices.
{
        local(res);
        res=w[1];
        for(i=2,length(w),res=res*w[i]);
        res
}
normk(a)=      \\ computes max{ |s(a)| }, s running over all 
               \\ embeddings s:Q[A]/(POL) --> C.
{
        local(kk,mx);
        kk=2*KLEVEL+4;mx=0;
        for(i=0,kk/2-1,if(gcd(1+2*i,kk) == 1, 
        mx=max(mx,abs(subst(lift(a),A,exp(2*Pi*(2*i+1)*I/kk))))));
        mx
}
cpx(t)=    \\ computes the characteristic polynomial of t over
           \\ the field Q[A]/(POL) in Q[A]/(POL)[x].
{
        charpoly(t,x,2)    \\ with option 2: using Hessenberg form.
}
mult_sparse(t,s)=   \\ t tabular matrix, s sparse_matrix,
                    \\ the result is t*s in tabular format.
{
local(res,u,a,b,i);
         a=matsize(t)[1];b=matsize(t)[2];
         u=2;i=s[u][1];
         res=matrix(a,b);
         while(i > 0,
         res[,s[u][2]]=res[,s[u][2]]+s[u][3]*t[,i];u++;i=s[u][1]);
         res
}
s2t(s)=      \\ restores the sparse matrix in tabular format.
{
	mult_sparse(matid(s[1][1]),s)
}
t2s(t)=      \\ transforms tabular to sparse format.
{               
        local(count,res,dimi,dimj);
        dimi=matsize(t)[1];dimj=matsize(t)[2];
        res=vector(dimi*dimj+2);
        res[1]=[dimi,dimj,"M"];
        count=1;
        for(i=1,dimi,
                for(j=1,dimj,
                        if( t[i,j] == 0,  , 
                                count++;res[count]=[i,j,t[i,j]];)));
        res=vector(count++,i,res[i]);
        res[count]=[0,0,0];
        res
}
slw(m)=    \\ a 2x2 integral, det=1,  matrix m is written as product of
           \\ a=Mat([1,1;0,1]);b=Mat([1,0;-1,1]);1/a;1/b.
{
        local(m1,a,b,g,tr,j,res,limitation);
        if( matdet(m) != 1, print("WARNING: Not in SL");return);
        a=Mat([1,1;0,1]);b=Mat([1,0;-1,1]);
        g=[a,1/a,b,1/b];
        m1=m;
        tr=trace(m1*mattranspose(m1));
	limitation=1000;      \\ not really a limitation
        res=vector(min(limitation,6+tr));   
        z=0;
        while(tr>2,
                trg=vector(4,i,trace(m1/g[i]*mattranspose(m1/g[i]) ) );
                for(i=1,4,if(trg[i]<tr,j=i;tr=trg[i]));
           z++;
	   if(z > limitation, 
	   return(print("WARNING: increase limitation in slw()")));
           res[z]=g[j];
           m1=m1/g[j];
           );

	   if(m1[1,1] == 0 && m1[1,2] == 1, z++;res[z]=a;
	                                    z++;res[z]=b;
					    z++;res[z]=a);
	   if(m1[1,1] == 0 && m1[1,2] == -1, z++;res[z]=1/a;
	                                     z++;res[z]=1/b;
					     z++;res[z]=1/a);
	   if(m1[1,1] == -1, z++;res[z]=a;
                             z++;res[z]=b;
                             z++;res[z]=a;
                             z++;res[z]=b;
                             z++;res[z]=a;
                             z++;res[z]=b);				    
	   vector(z,i,res[z+1-i]);
}
eval_sl(m,k,inp)=   \\ evaluates the TQFT so action of m \in SL(2,Z)
                        \\ on torus with one puncture, 
			\\ level k, input color inp.
			\\ if k>0 in "so" with k odd, else in "su".
{
        local(v,aa,bb,res);
        if( matdet(m) != 1, print("WARNING: Not in SL");return);
	v=slw(m);
        if(k>0,init_so(k,0),init_su(-k,0)); 
	init_boom([0],inp);
	if(matsize(BoomBasis)[2]==0,return(matrix(0,0)));
        aa=twA(1); bb=twB(1);
        res=matid(matsize(aa)[1]);
        for(i=1,length(v),if(v[i][1,2] == 1, res=res*aa);
                  if(v[i][1,2] == -1, res=res/aa);
		  if(v[i][2,1] == 1, res=res/bb);
		  if(v[i][2,1] == -1, res=res*bb);
		 );
		 res
}