version = "1.02";
(**

   fu.m            version 1.02  October 27, 2015
   
   Author: S. Orevkov   orevkov@math.ups-tlse.fr

   Computations for the paper [fu.pdf]
   "Markov trace on the Funar algebra", arXiv:1206.0765
   The output is in Macaulay2 format.

   This program is called either from a Mathematica session by the command

   <<fu.m

   or from UNIX command line by the command

   math <fu.m

   It creates a file  fu_dataL.m2  in the Macaulay2 format with the
   definitions of the matrix L0 (the generators of $\bar L^{0}$)
   and the matrices  rhoL_0,...,rhoL_8

   It creates also a file  fu_dataL.sing  in the Singular format
   with the definitions of the same matrices.

   Elements of the module $AF_3^+ \otimes AF_3^+$ are represented by
   linear combinations with coefficients in Z[a,b,u,v] of expressions of
   the form x[1**{i_1,..,i_k}**{j_1,..,j_m}**1] where i_1,.. and j_1,..
   are integers belonging to the set {-2,-1,1,2}. Such an expression
   represents  $X \otimes Y$  where

     $X = x_{|i_1|}^{sign(i_1)} \dots x_{|i_k|}^{sign(i_k)}$,
     $Y = x_{|j_1|}^{sign(j_1)} \dots x_{|j_m|}^{sign(j_m)}$.

   The list of rools `red' when applied to an element of
   $AF_3^+ \otimes AF_3^+$, computes the image of the element under
   the mapping $pr \circ (\tau_3 \otimes \tau_3)$ where
   pr  is the projection of the module $N = AF_2^+ \otimes AF_2^+$
   onto its submodule $\bar N$ generated by:

   e_0 = 1 \otimes 1,
   e_1 = x_1^{-1} \otimes x_1^{-1}
   e_2 = x_1^{-1} \otimes x_1
   e_3 = x_1      \otimes x_1^{-1}
   e_4 = x_1      \otimes x_1

   The projection  pr  maps other generators of  N  as follows:

   pr( x_1^{-1} \otimes 1 ) = pr( 1 \otimes x_1^{-1} ) = v e_0
   pr( x_1      \otimes 1 ) = pr( 1 \otimes x_1      ) = u e_0

**)

red={
x[X_**{ 2, 1, 2}**Y_]->u*x[X**{ 1, 1}**Y], (* { 1, 2, 1} -> u*{ 1, 1} *)
x[X_**{ 2, 1,-2}**Y_]->u*x[X**{-1, 1}**Y], (* {-1, 2, 1} -> u*{-1, 1} *)
x[X_**{-2, 1, 2}**Y_]->u*x[X**{ 1,-1}**Y], (* { 1, 2,-1} -> u*{ 1,-1} *)
x[X_**{ 2,-1,-2}**Y_]->v*x[X**{-1, 1}**Y], (* {-1,-2, 1} -> v*{-1, 1} *)
x[X_**{-2,-1, 2}**Y_]->v*x[X**{ 1,-1}**Y], (* { 1,-2,-1} -> v*{ 1,-1} *)
x[X_**{-2,-1,-2}**Y_]->v*x[X**{-1,-1}**Y], (* {-1,-2,-1} -> v*{-1,-1} *)
x[X_**{ 2,-1, 2}**Y_]->(
       (2a-b^2) x[X**{}**Y]
     -          (x[X**{ 1}**Y] + x[X**{ 2}**Y])
     - (a^2-b)  (x[X**{-1}**Y] + x[X**{-2}**Y])
     + b        (x[X**{ 1, 2}**Y] + x[X**{ 2, 1}**Y])
     + a        (x[X**{ 1,-2}**Y] + x[X**{ 2,-1}**Y] +
                 x[X**{-1, 2}**Y] + x[X**{-2, 1}**Y])
     + (a*b-1)  (x[X**{-1,-2}**Y] + x[X**{-2,-1}**Y])
     - a         u*x[X**{ 1, 1}**Y]
     -          (u*x[X**{ 1,-1}**Y] + u*x[X**{-1, 1}**Y] +
                 v*x[X**{ 1, 1}**Y])
     - b        (v*x[X**{-1, 1}**Y] + v*x[X**{ 1,-1}**Y])
     + (a-b^2)   v*x[X**{-1,-1}**Y]
),
x[X_**{-2, 1,-2}**Y_]->(
       (2b-a^2) x[X**{}**Y]
     -          (x[X**{-1}**Y] + x[X**{-2}**Y])
     - (b^2-a)  (x[X**{ 1}**Y] + x[X**{ 2}**Y])
     + a        (x[X**{-1,-2}**Y] + x[X**{-2,-1}**Y])
     + b        (x[X**{-1, 2}**Y] + x[X**{-2, 1}**Y] +
                 x[X**{ 1,-2}**Y] + x[X**{ 2,-1}**Y])
     + (a*b-1)  (x[X**{ 1, 2}**Y] + x[X**{ 2, 1}**Y])
     - b         v*x[X**{-1,-1}**Y]
     -          (v*x[X**{-1, 1}**Y] + v*x[X**{ 1,-1}**Y] +
                 u*x[X**{-1,-1}**Y])
     - a        (u*x[X**{ 1,-1}**Y] + u*x[X**{-1, 1}**Y])
     + (b-a^2)   u*x[X**{ 1, 1}**Y]
),
x[X_**{ 1, 1}**Y_]->( a*x[X**{ 1}**Y] - b*x[X**{}**Y] + x[X**{-1}**Y] ),
x[X_**{ 2, 2}**Y_]->( a*x[X**{ 2}**Y] - b*x[X**{}**Y] + x[X**{-2}**Y] ),
x[X_**{-1,-1}**Y_]->( b*x[X**{-1}**Y] - a*x[X**{}**Y] + x[X**{ 1}**Y] ),
x[X_**{-2,-2}**Y_]->( b*x[X**{-2}**Y] - a*x[X**{}**Y] + x[X**{ 2}**Y] ),
x[X_**{ 1,-1}**Y_]->x[X**{}**Y],
x[X_**{ 2,-2}**Y_]->x[X**{}**Y],
x[X_**{-1, 1}**Y_]->x[X**{}**Y],
x[X_**{-2, 2}**Y_]->x[X**{}**Y],
x[X_**{p_,q_,0 }**Y_]->x[X**{p,q}**Y],
x[X_**{p_,0 ,q_}**Y_]->x[X**{p,q}**Y],
x[X_**{0 ,p_,q_}**Y_]->x[X**{p,q}**Y],
x[X_**{0 ,p_}**Y_]->x[X**{p}**Y],
x[X_**{p_,0 }**Y_]->x[X**{p}**Y],
x[X_**{ 0}**Y_]->  x[X**{}**Y],
x[X_**{ 2}**Y_]->u*x[X**{}**Y],
x[X_**{-2}**Y_]->v*x[X**{}**Y],
x[X_**{ 1, 2}**Y_]->u*x[X**{ 1}**Y],
x[X_**{ 1,-2}**Y_]->v*x[X**{ 1}**Y],
x[X_**{-1, 2}**Y_]->u*x[X**{-1}**Y],
x[X_**{-1,-2}**Y_]->v*x[X**{-1}**Y],
x[X_**{ 2, 1}**Y_]->u*x[X**{ 1}**Y],
x[X_**{-2, 1}**Y_]->v*x[X**{ 1}**Y],
x[X_**{ 2,-1}**Y_]->u*x[X**{-1}**Y],
x[X_**{-2,-1}**Y_]->v*x[X**{-1}**Y],
x[1**{  }**{  }**1]->  e[0],
x[1**{  }**{ 1}**1]->u*e[0],
x[1**{ 1}**{  }**1]->u*e[0],
x[1**{  }**{-1}**1]->v*e[0],
x[1**{-1}**{  }**1]->v*e[0],
x[1**{-1}**{-1}**1]->  e[1],
x[1**{-1}**{ 1}**1]->  e[2],
x[1**{ 1}**{-1}**1]->  e[3],
x[1**{ 1}**{ 1}**1]->  e[4]}

(****  L0 is the list of generators of the module $\bar L^{(0)}$  ****)

L0=Union[Flatten[Table[
  Expand[(x[1**{2e1,e2,2e3}**{e4}**1]-x[1**{e2}**{2e3,e4,2e1}**1])//.red],
{e1,-1,1,2},{e2,-1,1},{e3,-1,1,2},{e4,-1,1}]]];
L0=Complement[L0,{0}];


(****   rhoL[[i]] is the list $\rho_i(e_0), \dots, \rho_i(e_4)$   ****)
(****   Each element of this list is represented as a polynomial  ****)
(****   in  a,b,u,v,e[0],..,e[4]  linear in  e[0],..,e[4]         ****)

eps={{0,0},{-1,-1},{-1,1},{1,-1},{1,1}};
de={{0,0},{0,-1},{0,1},{-1,0},{1,0},{-1,-1},{-1,1},{1,-1},{1,1}}
rhoL=Table[
  Table[
    Expand[ x[1**{de[[i,1]]}**{2eps[[j,1]],de[[i,2]],2eps[[j,2]]}**1]//.red ],
  {j,5}],
{i,9}];


(**** Check the invariance under the conjugation ****)

conj=Function[A,
  (A//.{a->aa,b->bb,u->uu,v->vv,e->ee})//.
    {aa->b,bb->a,uu->v,vv->u,ee[0]->e[0],
      ee[1]->e[4],ee[2]->e[3],ee[3]->e[2],ee[4]->e[1]}
];
T={{1,0,0,0,0},
   {0,0,0,0,1},
   {0,0,0,1,0},
   {0,0,1,0,0},
   {0,1,0,0,0}};
Sort[Table[ Expand[conj[rhoL[[i]]].T], {i,9} ]] - Sort[rhoL]
Sort[conj[L0]] - Sort[L0]


(**** Output to  fuL_data.m2 (for Macaulay2) ****)

mat={
  e[0]->{1,0,0,0,0},
  e[1]->{0,1,0,0,0},
  e[2]->{0,0,1,0,0},
  e[3]->{0,0,0,1,0},
  e[4]->{0,0,0,0,1}
};
L0m2=Expand[L0//.mat];
rhoLm2=Expand[rhoL//.mat];

out = OpenWrite["fu_dataL.m2",PageWidth->75];
WriteString[out,
"-- This file is created by the Mathematica program  fu.m  version ",version,"\n",
"-- called from a Mathematica session by the command\n--\n-- << fu.m\n--\n",
"-- Computations for the paper arXiv:1206.0765\n--\n",
"-- This file contains the definitions of:\n",
"-- L_0  is the matrix of generators of the module $\\bar L^{(0)}$\n",
"-- rhoL_i, i=0,..,8,  are the matrices of the operators $\\rho_\\delta$\n--\n",
"L_0=matrix(\n"];

Transpose[L0m2] >>> "fu_dataL.m2"
WriteString[out,")\n"];

Do[
  WriteString[out,"\nrhoL_",i-1,"=matrix(\n"];
  Transpose[rhoLm2[[i]]] >>> "fu_dataL.m2";
  WriteString[out,")\n"],
{i,9}];
Close[out];


(**** Output to  fuL_data.sing (for Singular) ****)

eZ={e[0]->Z[1],e[1]->Z[2],e[2]->Z[3],e[3]->Z[4],e[4]->Z[5]};
L0sing=Expand[L0//.eZ];
rhoLsing=Expand[rhoL//.eZ];

out = OpenWrite["fu_dataL.sing",PageWidth->75];
WriteString[out,
"// This file is created by the Mathematica program  fu.m  version ",version,"\n",
"// called from a Mathematica session by the command\n//\n// << fu.m\n//\n",
"// Computations for the paper arXiv:1206.0765\n//\n",
"// This file contains the definitions of:\n",
"// L0  is the matrix of generators of the module $\\bar L^{(0)}$\n",
"// rhoL(i), i=0,..,8,  are the matrices of the operators $\\rho_\\delta$\n//\n",
"module L0=\n"];

Do[
  WriteString[out,
    ToString[L0sing[[i]],FormatType->InputForm],
    If[i<Length[L0],",\n",";\n\n"]
  ],
{i,Length[L0]}]

Do[
  WriteString[out,"module rhoL(",ToString[i-1],")=\n"];
  Do[
    WriteString[out,
      ToString[rhoLsing[[i,j]],FormatType->InputForm],
      If[j<5,",\n",";\n\n"]
    ],
  {j,5}],
{i,9}]
Close[out];
