(** mtfa.m                                          version 0.14 **)
(** Last update: January 14, 2016                                **)
(**                                                              **)
(** Computation of Markov Trace on Funar Algebra                 **)
(** and the associated link invariants                           **)
(**                                                              **)
(** Author: Stepan Orevkov  orevkov@math.ups-tlse.fr             **)
(**                                                              **)
(** This file contains definitions of the following functions    **)
(** whose meaning is explained below:                            **)
(**                                                              **)
(** mtfa, mtfaP, mtfaPQ, mtfaPZ, mtfaPZlexico, mtfaP37,          **)
(** mtfaP37u2, GBQreduc, GBZreduc                                **)
(**                                                              **)
(** Input braids br are in the format of the package KnotTheory  **)
(** Example: BR[3,{2,-1,2,-1}]  is a braid representing the      **)
(**                             the figure-eight knot            **)
(**                                                              **)
(** The names of the variables  a, b, u, v, U, V  are fixed;     **)
(** U  and  V  are used instead of  $\bar u$  and  $\bar v$.     **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** mtfa[br,a,b,u,v]    Markov Trace on the Funar Algebra        **)
(**                                                              **)
(** input:  br          braid (in the format as explained above) **)
(**         a,b,u,v     the parameters (any expressions)         **)
(**                                                              **)
(** output: t(br) as a polynomial in  a,b,u,v                    **)
(**         the result is not yet a link invariant               **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** mtfaP[br,a,b,u,v,U,V]  raw P(br) where:                      **)
(**                                                              **)
(**           P(br) = u^((-n-e+1)/2)*v^((-n+e+1)/2)*t(br)        **)
(**                   when the number of link components is odd  **)
(**                   (in particular, for knots)                 **)
(**           P(br) = u^((-n-e)/2)*v^((-n+e)/2)*t(br)            **)
(**                   when the number of link components is even **)
(**                                                              **)
(** input:  br          braid (in the format as explained above) **)
(**         a,b,u,v     the parameters (any expressions)         **)
(**         U,V         stand for 1/u and 1/v respectively       **)
(**                                                              **)
(** output: P(br) as a polynomial in  a,b,u,v,U,V                **)
(**         the result is not yet a link invariant as it is but  **)
(**         it will be a link invariant after a ruduction mod I  **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** mtfaPQ[br]  normalized MTFA reduced mod I with b=0 and with  **)
(**             coefficients in Q (the same as in F_p where      **)
(**             a prime  p  is not in {2,3,7,37})                **)
(**                                                              **)
(** input:  either a braid  br  (the format is explained above)  **)
(**         or a precomputed Markov trace  t                     **)
(**                                                              **)
(** output: If the argument is a braid  br, then the output is   **)
(**         P(br) reduced mod I as a polynomial in  a,u          **)
(**         the reduction is computed with respect to the        **)
(**         lexicographic order with  V > U > v > u > a          **)
(**         (the reduced Groebner base with respect to this      **)
(**         order is given in the paper). So, the result is      **)
(**         a linear combination of the monomials:               **)
(**                                                              **)
(**         1, a^3, a^6, a^2*u, a^5*u, a*u^2, a^4*u^2, u^3       **)
(**                                                              **)
(**         If the argument is not of the form BR[...], then it  **)
(**         is interpreted as a precomputed Markov trace  t.     **)
(**         This means that the result of  mtfaPQ[br]  is        **)
(**         equivalent to  mtfaPQ[ mtfaP[br,a,0,u,v,U,V] ]       **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** mtfaPZ[br]  normalized MTFA reduced mod I with a=b=0  with   **)
(**             coefficients in  Z                               **)
(**                                                              **)
(** input:  the same as in  mtfaPQ                               **)
(**                                                              **)
(** output:  Similar to  mtfaPQ  but the reduction is computed   **)   
(**          modulo  I(0,0;Z). The monomial order is             **)
(**          V > U > w > v' > u'  where  u = w u' and  v = w v'  **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** mtfaPZlexico[br]  the same as mtfaPZ but the monomial order  **)
(**                   is  V < U < v < u                          **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** mtfaP37[br] normalized MTFA reduced mod I with b=0 and with  **)
(**             coefficients in GF(37) = Z/37Z                   **)
(**                                                              **)
(** input:  the same as in  mtfaPQ                               **)
(**                                                              **)
(** output:  Similar to  mtfaPQ  but the reduction is computed   **)   
(**          modulo  I(a,0;Z/37Z[a]). So, the result is          **)
(**         a linear combination of the monomials:               **)
(**                                                              **)
(**         1, a^3,a^6, a^2*u,a^5*u, a*u^2,a^4*u^2, u^3, a^2*u^4 **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** mtfaP37u2[br] MTFA mod 37 at u^3=2                           **)
(**                                                              **)
(** input:  the same as in  mtfaPQ                               **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** GBQreduc[G,f] reduction of a polynomial  f  mod a Groebner   **)
(**               base G (a list of polynomials) over Q.         **)
(**               The term order is always reverse lexicographic **)
(**               where the variables are orderd in the usual    **)
(**               alphabetic order used in Mathematica:          **)
(**               a < A < b < B < ... < z < Z                    **)
(**                                                              **)
(**--------------------------------------------------------------**)
(** GBZreduc[G,f] reduction of a polynomial  f  mod a Groebner   **)
(**               base G (a list of polynomials) over Z.         **)
(**               The term order is always reverse lexicographic **)
(**               where the variables are orderd in the usual    **)
(**               alphabetic order a < A < b < B < ... < z < Z   **)
(**                                                              **)
(**--------------------------------------------------------------**)
mtfa=Function[{br,alpha,beta,u,v},
Module[{a,b,n=br[[1]],rel,B,B0,i,j,k,m,nm,i0,ii,jj,ij,len,aa,bb,xSave=x,ySave=y,
  yXy={{},{1},{2},{-1},{-2},{1,2},{2,1},{1,-2},{2,-1},{-1,2},{-2,1},{-1,-2},{-2,-1},
       {1,2,1},{-1,2,1},{1,-2,1},{1,2,-1},{-1,-2,1},{1,-2,-1},{-1,-2,-1}},
  YxY,coef,cyXy,cYxY,found,foundij,c,f,rhs,crhs,(*stack,*)
  mtfaDebug=False}, Clear[x,y];
(* *)
  YxY=-yXy;
  coef={2a-b^2,-1,-1,b-a^2,b-a^2,b,b,a,a,a,a,a*b-1,a*b-1,-a,-1,-1,-1,-b,-b,a-b^2};
  cyXy=Expand[coef/.{a->alpha,b->beta}]; cYxY=Expand[coef/.{a->beta,b->alpha}];
  Do[If[cyXy[[i]]===0,yXy=Delete[yXy,i];cyXy=Delete[cyXy,i]],{i,20,1,-1}];
  Do[If[cYxY[[i]]===0,YxY=Delete[YxY,i];cYxY=Delete[cYxY,i]],{i,20,1,-1}];
  If[mtfaDebug,Print["yXy=",yXy];Print["cyXy=",cyXy];Print["YxY=",YxY];Print["cYxY=",cYxY]];
  rel=Flatten[Table[{x_**i**(-i)**y_->x**y,x_**(-i)**i**y_->x**y},{i,n-1}]];
  If[alpha===0 && beta===0, rel=Join[rel,Flatten[
    Table[{x_**i**i**y_->x**(-i)**y,x_**(-i)**(-i)**y_->x**i**y},{i,n-1}]
  ]]];
  Do[j=i+1;
    rel=Union[rel,{
      x_**  j **  i **  j **y_->x**  i **  j **  i **y,
      x_**(-j)**  i **  j **y_->x**  i **  j **(-i)**y,
      x_**  j **  i **(-j)**y_->x**(-i)**  j **  i **y,
      x_**(-j)**(-i)**  j **y_->x**  i **(-j)**(-i)**y,
      x_**  j **(-i)**(-j)**y_->x**(-i)**(-j)**  i **y,
      x_**(-j)**(-i)**(-j)**y_->x**(-i)**(-j)**(-i)**y
    }],
  {i,n-2}];
  Do[Do[
    rel=Union[rel,{
      x_**  j **  i **y_->x**  i **  j **y,
      x_**  j **(-i)**y_->x**(-i)**  j **y,
      x_**(-j)**  i **y_->x**  i **(-j)**y,
      x_**(-j)**(-i)**y_->x**(-i)**(-j)**y
    }],
  {j,i+2,n-1}],{i,n-3}];
 (**=============================================================**
  **                                                             **
  **    We transform the braid into a noncommutative product,    **
  **    e.g. BR[3,{1,-2,1,-2}] becomes 0**1**(-2)**1**(-2)**0    **
  **                                                             **
  **=============================================================**)
  B=Join[{0},br[[2]],{0}]; B[[0]]=NonCommutativeMultiply;
  stack={{1,B}}; f=0; mtfaMaxStackLength=1; mtfaCount=0;
  If[mtfaDebug,br>>mtfa.dbg];
  While[(len=Length[stack])>0,
    If[len>mtfaMaxStackLength, mtfaMaxStackLenth=len]; mtfaCount++;
    If[mtfaDebug,Print["======="];Print["f=",f];Print["stack=",stack]];
    If[mtfaDebug,"======">>>mtfa.dbg; f>>>mtfa.dbg; stack>>>mtfa.dbg];
    c=stack[[-1,1]]; B=stack[[-1,2]]//.rel;
    If[mtfaDebug,Print[" B=",B]];
    If[mtfaDebug," B=">>>mtfa.dbg; B>>>mtfa.dbg];
    If[ (len=Length[B])==2, f=Expand[f+c]; stack=Delete[stack,-1]; Continue[] ];
   (**================================**
    **                                **
    **   Step 1. Try to destabilize   **
    **                                **
    **================================**)
    m=0; i0=0; nm=0;
    Do[If[(k=Abs[B[[i]]])>m, m=k;nm=1;i0=i, If[k==m,nm++]],{i,2,len-1}];
    If[nm==1,stack[[-1]]={Expand[If[B[[i0]]>0,u,v]c],Delete[B,i0]};Continue[]];
   (**===========================================**
    **                                           **
    **   Step 2. Try to apply a braid relation   **
    **                                           **
    **===========================================**)
    found=False; foundij=False;
    Do[
      If[Abs[B[[i]]]-(k=Abs[B[[i+1]]])==1,
        Do[
          If[Abs[B[[j]]]==k || Abs[B[[j]]]==k+2, Break[] ];
          If[Abs[B[[j]]]==k+1,
            If[ B[[i]]==B[[j]],
              If[ (sgn=Sign[B[[i]]])==Sign[B[[i+1]]], found=True;
                If[mtfaDebug,Print["case yxy or YXY"]];
                B[[i]]=k*sgn; B[[i+1]]=(k+1)sgn;
                B=Insert[Delete[B,j],k*sgn,i+2]; stack[[-1,2]]=B,
              (*Else*)
                foundij=True;ij={i,j}       (** to use in Step 4 **)
              ],
            (*Else*)
              If[mtfaDebug,Print["case yxY, Yxy, yXY, or YXy"]]; found=True;
              B[[i+1]]=Sign[B[[i+1]]](k+1);
              B[[i]]=Sign[B[[j]]]k;
              B=Insert[Delete[B,j],-B[[i]],i+2]; stack[[-1,2]]=B
            ];
            Break[]
          ],
        {j,i+2,len-1}];
        If[found,Break[]]
      ],
    {i,2,len-3}];
    If[found,Continue[]];
   (**===========================================**
    **                                           **
    **   Step 3. Try to apply xx -> ax - b + X   **
    **                     or XX -> bX - a + x   **
    **                                           **
    **===========================================**)
    If[Not[alpha===0 && beta===0], (** when a=b=0, this is in `rel' **)
      Do[
        If[ (k=B[[i]])==B[[i+1]], found=True;
          B[[i]]=-k; B=Delete[B,i+1]; stack[[-1,2]]=B;
          If[ k>0, aa=alpha;bb=beta, aa=beta;bb=alpha ];
          If[Not[aa===0], B[[i]]=k;     AppendTo[stack,{Expand[ c*aa],B}]];
          If[Not[bb===0], B=Delete[B,i];AppendTo[stack,{Expand[-c*bb],B}]];
          Break[]
        ],
      {i,2,len-2}];
      If[ found, Continue[] ];
    ];
   (**====================================================**
    **                                                    **
    **   Step 4. Try to apply the long Funar's relation   **
    **                                                    **
    **====================================================**)
    If[ Not[foundij], Return[ERROR] ];
    i=ij[[1]]; j=ij[[2]]; k=Abs[B[[i+1]]]; sgn=Sign[B[[i+1]]];
    If[mtfaDebug,Print["case yXy or YxY"]];
    rhs = If[B[[i]]>0,yXy,YxY]/.{2->k+1,-2->-k-1,1->k,-1->-k};
    crhs = If[B[[i]]>0,cyXy,cYxY];
    stack=Delete[stack,-1];
    B=Delete[Delete[Delete[B,j],i+1],i];
    Do[B0=B;
      Do[B0=Insert[B0,rhs[[ii,jj]],i],{jj,Length[rhs[[ii]]]}];
      AppendTo[stack,{Expand[crhs[[ii]]c],B0}],
    {ii,Length[rhs]}];
  ]; (* EndWhile *) x=xSave; y=ySave;
  f
]];

mtfaP = Function[{br,a,b,u,v,U,V},Module[{t,n,e,pu,pv},
  t=mtfa[br,a,b,u,v];
  n=br[[1]]; e=Plus@@Sign[br[[2]]]; If[Mod[n+e,2]==0,n++];
  pu=(-n-e+1)/2; pv=(-n+e+1)/2;
  Expand[ If[pu>0,u^pu,U^(-pu)] * If[pv>0,v^pv,V^(-pv)] * t ]
]];

mtfaPQ = Function[br,Module[{t,G,g1,g2,g3,UU,VV,vv,aSave,uSave,vSave,USave,VSave},
  aSave=a;uSave=u;vSave=v;USave=U;VSave=V;Clear[a,u,v,U,V];
  If[Not[a===aSave],Print["WARNING: the variable `a' is redefined"]];
  If[Not[u===uSave],Print["WARNING: the variable `u' is redefined"]];
  g1=a^3+8; g2=2a^3+1; g3=3a^3+8;
  G=Expand[{g1*g2*g3,g1*g3(u-a),g3(6u^3-3a^2u+a^3+2),
    336u^4-792u^3*a+180u^2a^5+1272u^2a^2+846u*a^3+3264u-114a^7-1405a^4-3152a}];
  vv=-(336u^3a^2+216u^2a^3+2016u^2-432u*a^4-2112u*a-6a^8+53a^5+472a^2)/288;
  UU=-(170u^3a^5+528u^3a^2-18u^2a^6-80u^2a^3+384u^2-45u*a^7+56u*a^4+9a^8-30a^5+80a^2)/128;
  VV=-(-5880u^3a^4-5796u^2a^5-47712u^2a^2+14883u*a^6+110032u*a^3+118592u
       -8583a^7-78448a^4-130688a)/2688;
  t=If[br[[0]]===BR,mtfaP[br,a,0,u,v,U,V],br]/.{v->vv,U->UU,V->VV};
  v=vSave; U=USave; V=VSave; GBQreduc[G,t]
]];

mtfaPZlexico = Function[br,Module[{t,G,uSave,vSave,USave,VSave},
  uSave=u;vSave=v;USave=U;VSave=V;zUSave=zU;zVSave=zV;Clear[u,v,U,V,zU,zV];
  If[Not[u===uSave],Print["WARNING: the variable `u' is redefined"]];
  If[Not[v===vSave],Print["WARNING: the variable `v' is redefined"]];
  If[Not[U===USave],Print["WARNING: the variable `U' is redefined"]];
  If[Not[V===VSave],Print["WARNING: the variable `V' is redefined"]];
  G={U*u-1,V*v-1,v^3+v*u+u^3-3,V*U+5V*u^2+5U*v^2+5,V*u^3-3V+v^2+u,
     4u^3+4,4V-4u,4v+4u^2,4U+4u^2,16}/.{U->zU,v->zV};
  t=Expand[If[br[[0]]===BR,mtfaP[br,0,0,u,v,U,V],br]]/.{U->zU,v->zV};
  t=GBZreduc[G,t]/.{zU->U,zV->V}; zU=zUSave; zV=zVSave; t
]];

mtfaPZ = Function[br,Module[
  {t,n,e,G,su,uSave,vSave,wSave,USave,VSave,zUSave,zVSave,zWSave},
  uSave=u;vSave=v;wSave=w;USave=U;VSave=V;zUSave=zU;zVSave=zV;zWSave=zW;
  Clear[u,v,w,U,V,zU,zV,zW];
  If[Not[u===uSave],Print["WARNING: the variable `u' is redefined"]];
  If[Not[v===vSave],Print["WARNING: the variable `v' is redefined"]];
  If[Not[U===USave],Print["WARNING: the variable `U' is redefined"]];
  If[Not[V===VSave],Print["WARNING: the variable `V' is redefined"]];
  su={u->u*w,v->v*w,U->zU*zW,v->zV*zW};
  G={U*u-1,V*v-1,v^3+v*u+u^3-3,V*u^3-3V+v^2+u,V*U+5V*u^2+5U*v^2+5,
     4u^2+4v,4v*u-4,4v^2+4u,4U-4v,4V-4u,16};
  t=Expand[If[br[[0]]===BR,mtfaP[br,0,0,u,v,U,V],br]]/.su;
  t=GBZreduc[G,t]/.{zU->U,zV->V,w->1,zW->1};
  w=wSave; zU=zUSave; zV=zVSave; zW=zWsave; t
]];

mtfaP37 = Function[br,Module[
  {t,G,UU,VV,vv,aSave=a,uSave=u,vSave=v,USave=U,VSave=V},
  Clear[a,u,v,U,V];
  If[Not[a===aSave],Print["WARNING: the variable `a' is redefined"]];
  If[Not[u===uSave],Print["WARNING: the variable `u' is redefined"]];
  G={a^9+5a^6+2a^3-14,u*a^6-14u*a^3+9u-a^7+14a^4-9a,
     u^3*a^3+15u^3+18u*a^5+11u*a^2-6a^6+9a^3+5,
     u^5+2u^4*a+7u^3*a^2-9u^2*a^3+9u^2+6u*a^4+2u*a-12a^8+9a^5,37};
  vv=4u^4*a-15u^3*a^2+14u^2*a^3+16u^2+u*a^4+18u*a+4a^8-2a^5+2a^2;
  UU=3u^4*a+6u^3*a^2-6u^2*a^3+2u^2+10u*a^4-10u*a+2a^8-15a^5-14a^2;
  VV=-14u^4-16u^3*a+12u^2*a^5+11u^2*a^2-12u*a^3+u+12a^7+6a^4+a;
  t=Expand[If[br[[0]]===BR,mtfaP[br,a,0,u,v,U,V],br]];
  t=Expand[t/.{v->vv,U->UU,V->VV}];
  t=GBZreduc[G,t]; v=vSave;U=USave;V=VSave; t
]];

mtfaP37u2=Function[br,Module[{u3,aSave=a,uSave=u,vSave=v,USave=U,VSave=V,t},
  t=If[br[[0]]===BR,mtfaP[br,a,0,u,v,U,V],br];
  t=Mod[(PowerExpand[(t//.{v->-5u^2,a->-16u,U->-18u^2,V->11u})/.u->u3^(1/3)]/.u3->2),37];
  a=aSave;u=uSave;v=vSave;U=USave;V=VSave; t
]];


(****  Reduction of a polynomial ff using a Groebner base G over Q  ****)

GBQreduc=Function[{G,ff},
Module[{f=Expand[ff],f0,rc={True,0},i,reducDebug=False,momom,reducmonom},
 (* monom[f,i] -> The i-th monomial of f               *)
 (* When i=1, this is the leading monomial             *)
 (* Returns 0 if the numer of monomials is less than i *)
  monom=Function[{f,i},If[f[[0]]===Plus,If[i>Length[f],0,f[[-i]]],If[i>1,0,f]]];
  reducmonom=Function[{G,ff,ff0},
  Module[{f=ff,f0=ff0,cf0,denf0,i,g,lmg,lcg,deng,found=False},
    f0=(denf0=Denominator[f0])f0;
    cf0 = FactorList[f0][[1,1]];
    If[reducDebug, Print["f=",f] ];
    If[reducDebug, Print["f0,cf0,denf0=",{f0,cf0,denf0}] ];
    Do[g=G[[i]]; lmg=monom[g,1];
      If[PolynomialQ[Factor[f0/lmg]],
        lmg=(deng=Denominator[lmg])lmg;
        lcg=FactorList[lmg][[1,1]];
        f=Expand[f-g(f0/denf0)(deng/lmg)]; found=True; Break[]
      ],
    {i,Length[G]}];
    {found,f}
  ]];
  While[ rc[[1]] && Not[f===0], f=(rc=reducmonom[G,f,monom[f,1]])[[2]] ];
  If[ f[[0]] === Plus, i=1;
    While[ Not[(f0=monom[f,i])===0],
      f=(rc=reducmonom[G,f,f0])[[2]]; If[Not[rc[[1]]],i++]
    ]
  ];
  f
]];

(****  Reduction of a polynomial ff using a Groebner base G over Z  ****)

GBZreduc=Function[{G,ff},
Module[{f=Expand[ff],f0,rc={True,0},i,reducDebug=False,momom,reducmonom,count=0},
 (* monom[f,i] -> The i-th monomial of f               *)
 (* When i=1, this is the leading monomial             *)
 (* Returns 0 if the numer of monomials is less than i *)
  monom=Function[{f,i},If[f[[0]]===Plus,If[i>Length[f],0,f[[-i]]],If[i>1,0,f]]];
  reducmonom=Function[{G,ff,ff0},
  Module[{f=ff,f0=ff0,cf0,i,g,lmg,lcg,r,q,found=False},
    cf0 = FactorList[f0][[1,1]];
    If[reducDebug, Print["f=",f] ];
    If[reducDebug, Print["f0,cf0,found=",{f0,cf0,found}] ];
    Do[g=G[[i]]; lmg=monom[g,1];
      If[PolynomialQ[Factor[f0/lmg]],
        lcg=FactorList[lmg][[1,1]];
        If[(r=Mod[cf0,lcg])==cf0,Continue[]];
        f=Expand[f-g(f0/cf0)((cf0-r)/lmg)]; found=True; Break[]
      ],
    {i,Length[G]}];
    {found,f}
  ]];
  While[ rc[[1]] && Not[f===0], f=(rc=reducmonom[G,f,monom[f,1]])[[2]] ];
  If[ f[[0]] === Plus, i=1;
    While[ Not[(f0=monom[f,i])===0],
      f=(rc=reducmonom[G,f,f0])[[2]]; If[Not[rc[[1]]],i++]
    ]
  ];
  f
]];









  
