#include <stdio.h>
#include <stdlib.h>
char               *version = "1.02";
/**
   fu.c            Last update:      October 27, 2015

   Computations for the paper [fu.pdf]
   "Markov trace on the Funar algebra", arXiv:1206.0765
   The output is in Macaulay2 (default) or Singular (-S key) format.
   
   Author: S. Orevkov   orevkov@math.ups-tlse.fr

   We copmute the generators of $J_4^{(0)}$ (the matrix J4_0) and
   the 315x315 matrices rho_0,...,rho_30 of the linear
   operators $\rho_X \in End_A(AF_4^{red})$ for $X \in S_4$
   (see the details in fu.pdf)

   Compilation:

      gcc fu.c -o fu

   Usage:

      ./fu [-0] [-S] [>outputfile}

   The files "fu_data.m2" and "fu_data0.m2" are created by the commands

      ./fu >fu_data.m2
      ./fu -0 >fu_data0.m2

   The files "fu_data.sing" and "fu_data0.sing" are created by the commands

      ./fu -S >fu_data.sing
      ./fu -S -0 >fu_data0.sing
   
   The option -0 is used to set beta=0
   The option -S is used to switch the output format to Singular
   If ">outputfile" is not specified, then the stdout is used

   ---------------------------------------------------------------
   If the output is in the format of Macaulay2, then it has the
   following form:
   
   J_40=  <a matrix with 315 rows and 25 columns>
   rho_0=  <a matrix 315 x 315>
   rho_1=  <a matrix 315 x 315>
     . . . . . .
   rho_30= <a matrix 315 x 315>

   Each matrix is presented as a list of columns in the form
   
   col0 | col1 | ... ;
   
   and each column is represented as a linear combination of the
   base columns Z_{0}, Z_{1}, ... (Z_{i} is a column-vector with 1 at
   the i-th position and with zeros elsewhere).
   
   Before inputing these matrices in Macaulay2, the base columns Z_{i}
   shuld be defined e.g. like this:
   
      R = ZZ[v,u,a]; Z = map(R^315,R^315,1); load "fu_data0.m2";
   
   or

      R = ZZ[v,u,a,b]; Z = map(R^315,R^315,1); load "fu_data.m2";
   
   The mathematical meaning of Z_{i} is the i-th element of the base
   of the module $AF_4{red}$. This base is defined in the beginning
   part of main()

   By a technical reason, the matrix J_40 terminates by a zero column
   (since this is a matrix of generators of a module, this column does
   not play any role).
 
   If _sh is defined by the preprocessor command #define _sh, then the
   matrix of $pr_5\circ sh:AF_4^{red}\to A_4^{red}$ is also computed.
   Here $pr_5:AF_5^+\to AF_4^+$ is a projection whose kernel
   is spanned by $F_5^+ \setminus F_4^+$.

   ---------------------------------------------------------------
   If the output is in the format of Singular, then it has the
   following form:
   
   J40= <a matrix with 315 rows and 25 columns>
   rho(0)=  <a matrix 315 x 315>
   rho(1)=  <a matrix 315 x 315>
     . . . . . .
   rho[30]= <a matrix 315 x 315>

   Each matrix is presented as a list of columns in the form
   
   col1 , col2 , ... ;
   
   and each column is represented as a linear combination of the
   base columns Z[1], Z[2], ... (Z[i] is a column-vector with 1 at
   the i-th position and with zeros elsewhere).

   NB: rho(i) are numbered starting with 0 whereas
       Z[j]   are numbered starting with 1

   Before inputing these matrices in Singular, the base columns Z[i]
   shuld be defined e.g. like this:

      ring R=0,(v,u,a),dp; int b=0;
      matrix Z[315][315]; Z=Z+1;
      execute(read("fu_data0.sing"));
   
   or

      ring R=0,(v,u,b,a),dp;
      matrix Z[315][315]; Z=Z+1;
      execute(read("fu_data.sing"));

   The mathematical meaning of Z[i] is the same as in the case when
   the Macaulay2 format is used.

**/

/**======================== SECTION 0 ======================**/

#define YES 1
#define NO  0
#define NEW 1
#define ADD 0
#define    stack_size 1000000
#define ch_stack_size 1000000
#define clen 10
#define mdeg 11
#define mmdeg (mdeg*mdeg)

#define debug_     /* #define debug => debug on   */

#define _J40       /* #define _J40 => compute m0  */
#define _rho       /* #define _rho => compute rho */
#define _sh_       /* #define _sh  => compute sh  */

/** At the first stage, the elements of A=Z[u,v,a,b]   **/
/** are computed as polynomials in a,b,c,d,e,f,g,u,v   **/
/** where c,...,g are defined as in the comments below **/
/** each monomial in these variables is stored in the  **/
/** memory as 10 consecutive words (type int).         **/
/** the offset of the coefficient is 0 and the offset  **/
/** of the exponents of the variables are as follows:  **/

#define offs_a  1
#define offs_b  2
#define offs_c  3  /** c = 2a - b^2 **/
#define offs_d  4  /** d = 2b - a^2 **/
#define offs_e  5  /** e =  a - b^2 **/
#define offs_f  6  /** f =  b - a^2 **/
#define offs_g  7  /** g = ab - 1   **/
#define offs_u  8
#define offs_v  9

/** When a monomial in a,b,...,u,v is stored in the    **/
/** stack[], the exponent of a variable x (x=a,...,v)  **/
/** is stored at the offset offs_x w.r.t the stack     **/
/** pointer sp (see the comment before reduce() )      **/


/**======================== SECTION 1 ======================**/

int opt_0 = NO;    /**  opt_0 == YES <=> set beta = 0 **/
int opt_Sing = NO; /**  opt_Sing == YES <=> Singular  **/
int opt_H = NO;    /**  opt_H == YES <=> homogenize   **/

             /** Sintax features of the output format **/

	     /**                   Singular   M2      **/
char *cmt;   /** comment prefix:    "//"      "--"    **/
char *lb;    /**                    "["       "_{"    **/
char *rb;    /**                    "]"       "}"     **/
char *lp;    /**                    "("       "_"     **/
char *rp;    /**                    ")"       ""      **/
char *delim; /** columns delimiter  ","       "|"     **/
char *cmd;   /** command            "module " ""      **/

char *var = "abcdefguv";

int  *binom[ mdeg ]; /** binomial coefficients, i.e., **/
                     /** binom[n][k] = n!/(k!(n-k)!)  **/

int *sbinom[ mdeg ]; /** signed binomial coefficients **/
   /** i.e., sbinom[m][k] is the coefficient of X^k   **/
   /** in the expansion of (X-1)^m                    **/
   /** Thus  sbinom[m][k] = (-1)^(m-k)*binom[m][k]    **/

int *ce_coef[ mdeg ][ mdeg ];
   /**  ce_coef[p][q][ k ]  is the coefficient of X^k **/
   /**  in the expansion of (2-X)^p (1-X)^q           **/

int conjug[clen+1]={0, 2,1, 4,3, 6,5, 7, 9,8, 10};
   /**  the involution of   Z[a,b,c,d,e,f,g,u,v]      **/
   /**  which exchanges: a<->b, c<->d, e<->f, u<->v   **/

typedef unsigned char byte;
typedef unsigned int  word;

int     stack[   stack_size]; int    sp=0;
char ch_stack[ch_stack_size]; int ch_sp=0;


/**======================== SECTION 2 ======================**/

/**  S1,S2,S3,S4 are $S_1,...,S_4$ defined in  fu.pdf       **/
/**  A string (in the C-language sense) composed of the     **/
/**  letters a,b,c,d,A,B,C,D encodes a word in $x_i^{\pm1}$ **/
/**                                                         **/
/**  a,b,c,d are $x_1, ... , x_4$ and                       **/
/**  A,B,C,D are $x_1^{-1},...,x_4^{-1}$                    **/

byte *S1[3] ={"","a","A"};
byte *S2[7] ={"","b","ba","bA","B","Ba","BA"};
byte *S3[15]={"","c","cb","cba","cbA","cB","cBa","cBA",
                 "C","Cb","Cba","CbA","CB","CBa","CBA"};
byte *S4[31]={"","d","dc","dcb","dcba","dcbA","dcB","dcBa","dcBA",
                     "dC","dCb","dCba","dCbA","dCB","dCBa","dCBA",
                 "D","Dc","Dcb","Dcba","DcbA","DcB","DcBa","DcBA",
                     "DC","DCb","DCba","DCbA","DCB","DCBa","DCBA"};

/** lenS1[i] is the length of the string S1[i] **/
/** lenS2[i] is the length of the string S2[i] **/
/** lenS3[i] is the length of the string S3[i] **/
/** lenS4[i] is the length of the string S4[i] **/

int lenS1[3] ={0,1,1};
int lenS2[7] ={0,1,2,2,1,2,2};
int lenS3[15]={0,1,2,3,3,2,3,3,
                 1,2,3,3,2,3,3};
int lenS4[31]={0,1,2,3,4,4,3,4,4,
                   2,3,4,4,3,4,4,
                 1,2,3,4,4,3,4,4,
                   2,3,4,4,3,4,4};

/** The arrays R[], R_len[], R_coef[] and R_offs[] encode the **/
/** the right hand side of Funar relations -- relations (2)   **/
/** and (3) in the paper [fu.pdf]                             **/
/**                                                           **/
/** The values of R_offs[] are used in push()                 **/
/** the value R_offs[i]=10 is used there as a "dummy offset": **/
/** when it is used in push(), then it causes a modification  **/
/** of the array entry stack[sp+10] which is never used later **/


byte *R[19] = {"x","y", "X","Y", "xy","yx", "xY","yX","Xy","Yx", "XY","YX",
                          "xyx",  "Xyx","xYx","xyX",  "XYx","xYX",  "XYX"};
int R_len [19]={ 1, 1, 1,1, 2,2, 2,2,2,2, 2,2,  3,  3, 3, 3,   3, 3,  3};
int R_coef[19]={-1,-1, 1,1, 1,1, 1,1,1,1, 1,1, -1, -1,-1,-1,  -1,-1,  1};
int R_offs[19]={10,10, 6,6, 2,2, 1,1,1,1, 7,7,  1, 10,10,10,   2, 2,  5};
            /**  1  1  f f  b b  a a a a  g g   a   1  1  1    b  b   e **/


/** the array enc[] serves to reencode the words S1[i]...S4[i] **/
/** according the rules 'a'->1,...,'d'->4, 'A'->-1,...,'D'->-4 **/

int enc[256];

/**======================== SECTION 3 ======================**/

/** col is the array used in the computation of the matrices $\rho_X$ **/
/** and $J_4^{(0)}$. It represents a column of these matrices         **/
/** col[i] is the i-th entry of the column - an element of Z[a,b,u,v] **/
/** It is known a priori that each such entry is of the form          **/
/** A0 + A1*u + A2*v, so col[i][0],col[i]|1],col[i][2] are A0,A1,A2   **/
/** Each of Aj (j=0,1,2) is a polynomial in a and b                   **/
/**                                                                   **/
/** The array  col  represents an element of $AF_4^{red}$ of the form **/
/** Sum( col[j]*base4[j] )                                            **/
/**                                                                   **/
/** The polynomials  A0, A1, A2  are stored as polynomials in the     **/
/** variables C = b^2/a, D = a^2/b multiplied by monomials in a and b **/
/**                                                                   **/
/** This means that each of  Aj (j=0,1,2) is represented as           **/
/**                                                                   **/
/**  Aj = a^pow_a[i][j] * b^pow_b[i][j] * Sum_{k=0..deg_c[i][j]}      **/
/**              Sum_{l=0..deg_d[i][j]} col[i][j][k + mdeg*l]*C^k*D^l **/

int col[315][3][mdeg*mdeg];
int pow_a[315][3];
int pow_b[315][3];
int deg_c[315][3];
int deg_d[315][3];
int wk[mdeg*mdeg];

/** base4[i] is the i-th element of the base of $AF_4^{red}$ **/

byte heap[6*315];
byte *base4[315];
int base4_len[315];

char Y[10]; int lenY;


void err(char *s){ puts(s); exit(9); }
void nomem(void){ puts("no memory"); exit(9); }
void errdeg(void){err("degree");}

void col_init(void);
void reduc( int sp, int ch_sp, int n, int mode );
void init_binom( void );
void col_add( int sp, char *p, int n );
void push(int sp,int c,int coef,int offs);
int  find( char *p, int n );
int  col_print( void );

long long iter=0;

/**=========================================================**/

void print_stack(void){ int i;
  for( i=0; i<10; i++ )printf(" %3d",stack[i]);
  putchar('\n');
  for( i=0; i<10; i++ )printf(" %3d",stack[i+10]);
  putchar('\n');
  for( i=0; i<10; i++ )printf(" %3d",stack[i+20]);
  putchar('\n');
}

/**=========================================================**/

main( int argc, char *argv[] ){

  int i1,i2,i3,i4,i5,i6,n,i,j,k,m,count,zero;

  char *p,*q;

                /** Get arguments and print a header to the output **/
  for( i=1; i<argc; i++ ){
    if( argv[i][0]=='-' ){
      switch( argv[i][1] ){
        case '0': opt_0 = YES; break;
	case 'S':
	case 's': opt_Sing = YES; break;
	case 'H':
	case 'h': opt_H = YES; break;
      }
    }
  }

  cmt   = ( opt_Sing ? "//" : "--" );
  lb    = ( opt_Sing ? "["  : "_{" );
  rb    = ( opt_Sing ? "]"  : "}"  );
  lp    = ( opt_Sing ? "("  : "_"  );
  rp    = ( opt_Sing ? ")"  : ""   );
  delim = ( opt_Sing ? ","  : "|"  );
  cmd = (opt_Sing ? "module " : "" );

  printf("%s This file is created by the C-program fu.c version %s\n",cmt,version);
  printf("%s called by the command line\n%s\n%s  %s",cmt,cmt,cmt,argv[0]);
  if( argc>1 )printf(" %s",argv[1]);
  if( argc>2 )printf(" %s",argv[2]);
  if( argc>3 )printf(" %s",argv[3]);
  printf("\n%s\n%s Computations for the paper arXiv:1206.0765\n",cmt,cmt);
  printf("%s\n%s This file contains the definitions of:\n",cmt,cmt);

#ifdef _J40
  printf("%s  J4%s0  is the matrix of generators of the module $J_4^{(0)}$.\n",cmt,(opt_Sing?"":"_"));
#endif
#ifdef _rho
  printf("%s  rho%si%s, i=0..30  are the matrices of $\\rho_X$, $X\\in S_4$,\n",
          cmt,   lp,rp);
  printf("%s  rho%s0%s, corresponds to $X = 1$.\n",
          cmt,   lp,rp);
#endif
#ifdef _sh
  printf("%s  sh  is the matrix of $sh\\in End_A(AF_4^{red})$\n",cmt);
#endif
                          /** Initialization **/
  init_binom();

  enc[(word)'a']=1;  enc[(word)'A']=-1;
  enc[(word)'b']=2;  enc[(word)'B']=-2;
  enc[(word)'c']=3;  enc[(word)'C']=-3;
  enc[(word)'d']=4;  enc[(word)'D']=-4;

                           /**  Compute the base of $AF_4^{red}$ **/
  count=0; q=heap;
  for( i3=0; i3<15; i3++ ){
  for( i2=0; i2< 7; i2++ ){
  for( i1=0; i1< 3; i1++ ){
    p = ch_stack;
    base4[count] = q;
    for( i=0; i<lenS1[i1]; i++ ) *p++ = enc[*q++=S1[i1][i]];
    for( i=0; i<lenS2[i2]; i++ ) *p++ = enc[*q++=S2[i2][i]];
    for( i=0; i<lenS3[i3]; i++ ) *p++ = enc[*q++=S3[i3][i]];
    n = base4_len[count] = lenS1[i1] + lenS2[i2] + lenS3[i3];
    if(find(ch_stack,n)!=count)err("ierr0");/** check if find() is correct **/
    count++;
  }}}
                        /** The base of $AF_4^{red}$ is computed  **/

#ifdef _J40             /** compute the generators of $J_4^{(0)}$ **/
  printf("\n%sJ4%s0=\n",cmd,(opt_Sing?"":"_")); count=0; zero=YES;
  for( i1=-1; i1<=1; i1+=2 ){
  for( i2=-1; i2<=1; i2++  ){
  for( i3=-1; i3<=1; i3+=2 ){
  for( i4=-1; i4<=1; i4++  ){
  for( i5=-1; i5<=1; i5++  ){
  for( i6=-1; i6<=1; i6+=2 ){

    /** loop over all words $X_3 X_2 X_1$ as in [fu.pdf; Sect.2.3, (J1)]    **/
    /** i.e. over all words of the form                                     **/
    /** $x_3^{i_1} x_2^{i_2} x_1^{i_3} x_3^{i_4} x_2^{i_5} x_3^{i_6}$       **/
    /** where $i_1,i_3,i_6 \in \{\pm1\}$ and $i_2,i_4,i_5\in\{ -1, 0, 1 \}$ **/
    
    col_init(); p=ch_stack;


    stack[0]= 1; n=m=0; for( j=1; j<10; j++ )stack[j]=0;
          p[n++]=3*i1; if(i2)p[n++]=2*i2; p[n++]=3*i3;
    if(i4)Y[m++]=1*i4; if(i5)Y[m++]=2*i5; Y[m++]=3*i6;
    lenY=m; reduc(0,0,n,1);

    stack[0]=-1; n=m=0; for( j=1; j<10; j++ )stack[j]=0;
    Y[m++]=3*i1; if(i2)Y[m++]=2*i2; if(i4)Y[m++]=1*i4;
    p[n++]=3*i3; if(i5)p[n++]=2*i5;       p[n++]=3*i6;
    lenY=m; reduc(0,0,n,2);

    if( !zero ) puts(delim);
    zero=col_print();

  }}}}}}
 /** add the zero column (something should follow "|") **/
  puts(opt_Sing ? "0;" : "0*Z_{0};");
#endif

#ifdef _rho
                         /** computation of rho[0]..rho[30] **/
  for( k=0; k<31; k++ ){
    printf("%srho%s%d%s=\n",cmd,lp,k,rp);
    for( i=0; i<315; i++ ){
      stack[0]=1; for( j=1; j<10; j++ )stack[j]=0;
      col_init();
      p = ch_stack; q = S4[k];  m = lenS4[k];
      for( j=0; j<m; j++ ) p[j] = enc[q[j]];
      p = ch_stack+m; q = base4[i];  n = base4_len[i];
      for( j=0; j<n; j++ ){
       i1 = enc[ q[j] ];
       p[j] = ( i1>0 ? i1+1 : i1-1 );
      }
      reduc(0,0,n+m,0); col_print(); if(i<314)puts(delim);else puts(";");
    }
  }
#endif

#ifdef _sh
  printf("%ssh=\n",cmd);
  for( i=0; i<315; i++ ){int zero=NO;
    p = ch_stack;  q = base4[i];  n = base4_len[i];
    for( j=0; j<n; j++ ){
      i1 = enc[ q[j] ];
      p[j] = ( i1>0 ? i1-1 : i1+1 );
      if( p[j]==0 ){ zero=YES; break; }
    }
    if(zero) printf("Z0");
    else{
      printf("Z_{%d}",find(p,n));
    }
    if(i<314)puts("|");else puts(";");
  }
#endif

}
/**==================================================**/
/** Initialize col[] to represent the zero element   **/
/** of $AF_4^{red}$ (recall that the array col[]     **/
/** serves for storing elements of $AF_4^{red}$)     **/

void col_init(void){ int m,w,j;
  for( m=0; m<315; m++ ){
    for( w=0; w<3; w++ ){
      for( j=0; j<mmdeg; j++ )col[m][w][j] = 0;
      pow_a[m][w] = pow_b[m][w] = deg_c[m][w] = deg_d[m][w] = 0; 
    }
  }
}

/**============================================================**/
/** reduc() is a recursive procedure.                          **/
/**                                                            **/
/** Mathematical meaning:                                      **/
/** Compute                                                    **/
/**   $\tau_5(qX)$    (if mode==0)                             **/
/**   or $q r(r(X)Y)$ (if mode==1)                             **/
/**   or $q r(Xr(Y))$ (if mode==2)                             **/
/** and add the result to the element of $AF_4^{red}$          **/
/** represented by the array col[].                            **/
/**                                                            **/
/** $q$ is a monomial in a,b,c,d,e,f,g,u,v with integer coeff  **/
/** If mode=0, then $X\in F_5^+$;                              **/
/** If mode>0, then $X,Y\in F_4^+$                             **/
/**                                                            **/
/** Arguments:                                                 **/
/**    sp : the offset in stack[] of the monomial $q$, i.e.    **/
/**         stack[sp]   : the integer coefficient,             **/
/**         stack[sp+1] : the exponent of a,                   **/
/**         stack[sp+2] : the exponent of b,                   **/
/**           . . . . .                                        **/
/**         stack[sp+7] : the exponent of g,                   **/
/**         stack[sp+8] : the exponent of u,                   **/
/**         stack[sp+9] : the exponent of v,                   **/
/** ch_sp : the offset in ch_stack[] of $X$,                   **/
/**     n : the word length of $X$, i.e.                       **/
/**         $X$ is stored in ch_stack[0},...,ch_stack[n-1]     **/
/**                                                            **/
/** the word $Y$ is stored in Y[0],...,Y[lenY]                 **/

void reduc( int sp, int ch_sp, int n, int mode ){
  int i,j,c,c1,c2,abs_c,abs_c1,abs_c2,change=YES,n4,i4,wk,found;
  char *p = ch_stack + ch_sp;

#ifdef debug
  if((iter++)>6000)err("iter");
  printf("reduc: mode=%d sp=%d n=%d b=",mode,sp,n);
  for( j=0; j<n; j++ )printf(" %d",p[j]);
  putchar('\n');
#endif

  if( opt_0 ){
    if( stack[sp+offs_b] )return;  /** if deg_b > 0 **/
  }

  while( change ){ change=NO; n4=0;

#ifdef debug
  printf("  while: ");
  for( j=0; j<n; j++ )printf(" %d",p[j]);
  putchar('\n');
#endif

    for( i=0; i<n; i++ ){ c=p[i]; abs_c = ( c>0 ? c : -c );

#ifdef debug
  printf("    for: i=%d c=%d  br=",i,c);
  for( j=0; j<n; j++ )printf(" %d",p[j]);
  putchar('\n');
#endif

      if( abs_c == 4 ){ n4++; i4=i; }
      if( i == n-1 )break;
      c1=p[i+1];
      if( c1 == -c ){                          /** xX **/
        for( j=i+2; j<n; j++ ) p[j-2] = p[j];
        change = YES; n-=2; break;
      }
      if( c1 ==  c ){                          /** xx **/
        for( j=i+2; j<n; j++ ) p[j-1] = p[j];
	change = YES; n--; p[i] = -c;                /** X **/

        if( ch_sp + n + n >= ch_stack_size )nomem();
	for( j=0; j<n; j++ ) p[j+n] = p[j];
	push(sp,c,1,offs_a); p[i+n] = c;
	reduc( sp+clen, ch_sp+n, n, mode);          /** ax **/
	
        for( j=0;   j<i; j++ ) p[j+n  ]=p[j];
	for( j=i+1; j<n; j++ ) p[j+n-1]=p[j];
	push(sp,c,-1,offs_b);
	reduc( sp+clen, ch_sp+n, n-1, mode);        /** -b **/

        break;
      }

      abs_c1 = (c1>0 ? c1 : -c1);
      if( abs_c1 < abs_c - 1 ){                   /** xy = yx **/
        p[i]=c1; p[i+1]=c; change=YES; break;
      }
      if( i==n-2 )continue;

      if( abs_c1 != abs_c - 1 )continue;

      for( found=NO,j=i+2; j<n; j++ ){ c2=p[j];
        if( (abs_c2 = (c2>0 ? c2 : -c2)) == abs_c ){ change=found=YES; break; }
        if( abs_c2 > abs_c - 2 )break;
      }

      if( ! found )continue;

      if( ch_sp + n + n >= ch_stack_size )nomem();
      for( ; j>i+2; j-- ) p[j] = p[j-1];

      if( c2 == -c ){
        if( c>0 ){
	  if( c1==c-1 ){ p[i]=-c1; p[i+1]=c;  p[i+2]= c1; } /** yxY = Xyx */
	  else         { p[i]= c1; p[i+1]=c2; p[i+2]=-c1; } /** yXY = XYx */
	}
	else{ /** c<0 **/
	  if( c1==c2-1 ){p[i]= c1; p[i+1]=c2; p[i+2]=-c1;}  /** Yxy = xyX */
	  else          {p[i]=-c1; p[i+1]=c;  p[i+2]= c1;}  /** YXy = xYX */
	}
      }
      else{
        if( (int)c*c1 > 0 ){p[i]=p[i+2]=c1; p[i+1]=c;} /** yxy = xyx  YXY = XYX */
        else{                                       /** yXy = */
          int k,m; byte *q,*r;

          for( k=0; k<19; k++ ){ q=p+n; r=R[k]; m=R_len[k];
            enc[(word)'y']= c;  enc[(word)'Y']=-c;
            enc[(word)'x']=-c1; enc[(word)'X']=c1;
#ifdef debug
printf("      for k=%d y=%d (sp=%d)\n",k,enc[(word)'y'],sp);
#endif
            for( j=0;   j<i; j++ ) *q++ = p[j];
            for( j=0;   j<m; j++ ) *q++ = enc[ *r++ ];
            for( j=i+3; j<n; j++ ) *q++ = p[j];
            push(  sp, c, R_coef[k], R_offs[k] );
            reduc( sp+clen, ch_sp+n, n+m-3, mode );
          }

          for( j=i+3; j<n; j++ ) p[j-3] = p[j];
          stack[ sp + (c>0 ? offs_c : offs_d) ]++; n-=3;

          break;
        }
      }
    }   /** end of for( i=0; i<n; i++ ){   **/
    if( change )continue;
    if( n4 == 1 ){ change = YES;
      n--;  c = p[i4];
      for( i=i4; i<n; i++ ) p[i] = p[i+1];
      stack[ sp + (c>0 ? offs_u : offs_v) ]++;
    }
  }  /** end of while( change ){ **/
#ifdef debug
  print_stack();
#endif
  if( mode==0 ){ col_add( sp, p, n ); iter++; return; }
  if( mode==1 ){
    for( j=0; j<lenY; j++ )p[n+j]=Y[j];
    reduc( sp, ch_sp, n+lenY, 0 );
  }
  else{
    for( j=n-1; j>=0; j-- )p[j+lenY]=p[j];
    for( j=0; j<lenY; j++ )p[j]=Y[j];
    reduc( sp, ch_sp, n+lenY, 0 );
  }
}
/**======================= end of reduc() =====================**/

/**============================================================**/
/** col_add()                                                  **/
/**                                                            **/
/** Mathematical meaning: col += q*X where q is a monomial in  **/
/** a,b,...,u,v stored in stack[sp],...,stack[sp+9] (in the    **/
/** way as explained in the comment before reduc() )           **/
/** and X is an element of $F_4^{red}$ stored in p[1]...p[n-1] **/
/**                                                            **/
/** Here the monomial in a,b,...,u,v is expanded into a        **/
/** polynomial in C and D where C = b^2/a, D = a^2/b           **/
/** multiplied by the monomial a^(pow_a)*b^(pow_b).            **/
/** See the comments in SECTION 3 before main()                **/

void col_add( int sp, char *p, int n ){
  int i,j,k,m,w,xa,xb,xc,xd,xe,xf,xg,pa,pb,dc,dd,ya,yb,zc,zd,jj,
      corr_c,corr_d,corr_cd,c,z,
      *x,*col_mw,*col_mwj,*col_mwj_corr,*ce,*df,*g;

  x = stack+sp; c = *x;
  m = find(p,n); if(n != base4_len[m] )err("mmmm");

  xa = x[offs_a]; xb = x[offs_b]; xc = x[offs_c]; xd = x[offs_d];
  xe = x[offs_e]; xf = x[offs_f]; xg = x[offs_g];
  w = x[offs_u] + 2*x[offs_v];
  pa = pow_a[m][w]; pb = pow_b[m][w];
  dc = deg_c[m][w]; dd = deg_d[m][w];

#ifdef debug
  printf("col_add:: sp=%d m=%d w=%d  br =",sp,m,w);
  for( i=0; i<n; i++ ) printf(" %d",p[i]); putchar('\n');
  printf(
  "c=%d xa=%d xb=%d xc=%d xd=%d xe=%d xf=%d xg=%d  pa=%d pb=%d dc=%d dd=%d\n",
   c,   xa,   xb,   xc,   xd,   xe,   xf,   xg,    pa,   pb,   dc,   dd);
#endif

  /**  C = a^(-1)*b^2 **/
  /**  D = b^(-1)*a^2 **/

  col_mw = col[m][w];

  ya = xa+xc+xe; yb = xb+xd+xf;
  if( dc==0 ){ pow_a[m][w]=pa=ya; pow_b[m][w]=pb=yb; }

  if( (zc = 2*(yb-pb) + (ya-pa))%3 )err("mod3"); zc/=3;
  if( (zd = 2*(ya-pa) + (yb-pb))%3 )err("mod3"); zd/=3;

  corr_c = (zc < 0 ? -zc : 0);
  corr_d = (zd < 0 ? -zd : 0);
  corr_cd = corr_c + mdeg*corr_d;

#ifdef debug
printf("ya=%d yb=%d  zc=%d zd=%d  corr_c=%d corr_d=%d",
        ya,   yb,    zc,   zd,    corr_c,   corr_d);
#endif

  if( dc + corr_c > mdeg )err("mdeg");
  if( dd + corr_d > mdeg )err("mdeg");

#ifdef debug
  for( i=0; i<dc; i++ ){
    for( j=0; j<dd; j++ )printf(" %3d",col[m][w][i+mdeg*j]);
    putchar('\n');
  }
#endif

  if( corr_c || corr_d ){
    for( jj=(dd-1)*mdeg; jj>=0; jj-=mdeg ){
      col_mwj=col[m][w]+jj;
      col_mwj_corr = col_mwj + corr_cd;
      for( i=dc-1; i>=0; i-- ){
        col_mwj_corr[i] = col_mwj[i]; col_mwj[i] = 0;
      }
    }
  }

  dc += corr_c; dd += corr_d; zc += corr_c; zd += corr_d;

  if( (i = zc + xc + xe + xg + 1) > dc )dc=i;
  if( (i = zd + xd + xf + xg + 1) > dd )dd=i;
  if( dc > mdeg || dd > mdeg )err("mdeg");

  deg_c[m][w] = dc;  deg_d[m][w] = dd;

#ifdef debug
  printf("  dc=%d dd=%d\n",dc,dd);
#endif

  pow_a[m][w] -= (2*corr_d - corr_c);
  pow_b[m][w] -= (2*corr_c - corr_d);

#ifdef debug
  for( i=0; i<dc; i++ ){
    for( j=0; j<dd; j++ )printf(" %3d",col[m][w][i+mdeg*j]);
    putchar('\n');
  }
  putchar('\n');
#endif

  col_mw=col[m][w];
  g=sbinom[xg]; ce=ce_coef[xe][xc]; df=ce_coef[xf][xd];
  for( i=0; i<=xc+xe; i++ ){
    for( j=0; j<=xd+xf; j++ ){
      for( k=0; k<=xg; k++ ){
        z = c*ce[i]*df[j]*g[k];
        if( z > 0x3fffffff || z < -0x3fffffff )err("too big coef");
#ifdef degug
printf("(i,j)=(%d,%d) += %d\n", zc+i+k,zd+j+k,c*ce[i]*df[j]*g[k]);
#endif
        z = (col_mw[ zc+i+k + (zd+j+k)*mdeg ] += z);
        if( z > 0x3fffffff || z < -0x3fffffff )err("too big coef");
      }
    }
  }
#ifdef debug
  for( i=0; i<dc; i++ ){
    for( j=0; j<dd; j++ )printf(" %3d",col[m][w][i+mdeg*j]);
    putchar('\n');
  }

  printf("pa=%d pb=%d dc=%d dd=%d\n",
    pow_a[m][w],pow_b[m][w],deg_c[m][w],deg_d[m][w]);
  col_print(); putchar('\n');
#endif
}

/*===============================================*/
void init_binom( void ){ int i,j,k,m,p,q;

                          /** In the comments below "GF" means  **/
			  /**       "Generating Function"       **/
			  
                          /** GF: binom[m] = (1+X)^m            **/
                          /** GF:sbinom[m} = (1-X)^m            **/
                          /**      = sum binom[m][k] (-1)^k X^k **/
                          /** sbinom[m][k] = (-1)^k bimom[m][k] **/
  for( m=0; m<mdeg; m++ ){
    binom [m] = (int *)malloc( sizeof(int)*(m+1) );
    sbinom[m] = (int *)malloc( sizeof(int)*(m+1) );
    binom[m][0] = binom[m][m] = 1;
    for( k=1; k<m;  k++ ) binom[m][k] = binom[m-1][k-1] + binom[m-1][k];
    for( k=0; k<=m; k++ )sbinom[m][k] = (k&1 ? -1 : 1) * binom[m][k];
  }

                          /** GF:coef_ce[p][q] = (1-X)^p (2-X)^q **/
  for( p=0; p<mdeg; p++ ){for( q=0; q<mdeg; q++ ){
    ce_coef[p][q] = (int *)malloc( sizeof(int)*(p+q+1) );
    for( k=0; k<=p+q; k++ )ce_coef[p][q][k]=0;
    for( i=0; i<=p; i++ ){for( j=0; j<=q; j++ ){
      ce_coef[p][q][i+j] += ((sbinom[p][i]*sbinom[q][j])<<(q-j));
    }}
  }}

                       /** GF: sbinom[m} = (X-1)^m                **/
                       /**       = sum binom[m][k] X^k (-1)^(m-k) **/
                       /** sbimon[m][k] = old.sbonom[m][k] (-1)^m **/
  for( m=0; m<mdeg; m++ ){
    if( m&1 ){
      for( k=0; k<=m; k++ )sbinom[m][k] = -sbinom[m][k];
    }
  }
  
}

/**================================================**/
/** push qX to the stacks where q is a monomial in **/
/** a,b,...,u,v and X an element of $F_5^+$        **/
/** (see other comments above)                     **/

void push(int sp,int c,int coef,int offs){
  int i,*x,*y;
#ifdef debug
  printf("push sp=%d c=%d coef=%d offs=%d\n",sp,c,coef,offs);
#endif
  if( sp + clen + clen > stack_size )nomem();
  x = stack+sp; y = x+clen;
  for( i=0; i<clen; i++ ) *y++ = *x++;
  if( c < 0 )offs = conjug[offs];
  x[0] = coef*x[0]; x[offs]++;
}

/**=================================================**/
/** given a reduced 4-braid p[0]...p[n-1], find the **/
/** number j such that it is the j-th element of    **/
/** base4                                           **/

/** This subroutine may be treated as a black box   **/
/** because its correctness is checked at the       **/
/** initialization of base4[] in main()             **/

int find( char *p, int n ){ int i=0, res=0;

  if( n==0 )return 0;
  switch( p[0] ){
    case  1: res+=1; i=1; break;
    case -1: res+=2; i=1; break;
    default: i=0;
  }
  if( i == n )return res;
  switch( p[i] ){
    case -2: res+=9;
    case  2: res+=3; i++;
      if( i == n )return res;
      switch( p[i] ){
        case -1: res+=6; i++; break;
        case  1: res+=3; i++; break;
        default: break;
      }
    default: break;
  }
  if( i == n )return res;
  switch( p[i] ){
    case -3: res+=(7*21);
    case  3: res+=21; i++;
      if( i == n )return res;
      switch( p[i] ){
        case -2: res+=63;
        case  2: res+=21; i++;
          if( i == n )return res;
          switch( p[i] ){
            case -1: return res+42;
            case  1: return res+21;
            default: err("ierr1");
          }
        default: err("ierr2");
      }
    default: err("ierr3");
  }
}
/**====================================================**/
int len_uv[3]; int star_uv[3]; int c1_uv[3];

int col_print( void ){
  int first=YES,zero=YES,star,c,i,j,m,w,da,db,dc,dd,len,*col_mw;

  for( m=0; m<315; m++ ){ len=0;
    for( w=0; w<3; w++ ){
      dc=deg_c[m][w]; dd=deg_d[m][w]; len_uv[w]=0; col_mw=col[m][w];
      for( i=0; i<dc; i++ ){ 
        for( j=0; j<dd; j++ ){
          if( opt_0 ){
            if( pow_b[m][w] + i + i - j )continue;
          }
	  if( col_mw[i+mdeg*j] != 0 ){
	    len_uv[w]++;
	    if( w==0 )len++;
	  }
	}
      }
      if( w>0 && len_uv[w]>0 )len++;
    }
    if( len==0 )continue;
    zero = NO;

    if( len>1 ){if(!first)putchar('+'); putchar('('); first=YES;}
    for( w=0; w<3; w++ ){ dc=deg_c[m][w]; dd=deg_d[m][w]; col_mw=col[m][w];
      if( len_uv[w]==0 )continue;
      if( w>0 && len_uv[w]>1 ){
        if(!first)putchar('+'); putchar('('); first=YES;
      }
      for( i=0; i<dc; i++ ){ 
        for( j=0; j<dd; j++ ){
	  if( (c=col_mw[i+mdeg*j]) == 0 )continue;
	  da = pow_a[m][w] + j + j - i;
	  db = pow_b[m][w] + i + i - j;

          if( opt_0 && db )continue;

          star=NO;
	  if( (!first) && c>0 )putchar('+');
          if( c==-1 )putchar('-');
          else if( c!=1){ star=YES; printf("%d",c); }
          if(da>0){if(star)putchar('*');star=YES;putchar('a');if(da>1)printf("^%d",da);}
          if(db>0){if(star)putchar('*');star=YES;putchar('b');if(db>1)printf("^%d",db);}
	  if( !star ){
            if( (w==0 && len>1)||(w>0 && len_uv[w]>1) )putchar('1');
	  }
	  first=NO;
	}
      }
      if( w>0 ){
        if( len_uv[w]>1 ){star=YES; putchar(')');}
        if(star)putchar('*'); star=YES; putchar( w==1 ? 'u' : 'v' );
	first = NO;
      }
    }
    if( len>1 ){star=YES; putchar(')');}
    if(star)putchar('*');
    printf("Z%s%d%s",lb,(opt_Sing ? m+1 : m),rb); first=NO;
  }
  return zero;
}

/** Mathematica code to check that the entries of ce_coef{] do not exceed 30 bits

Do[
  Print[{n,N[(Max@@Table[Binomial[n,p]2^p,{p,n-1}])^2/2^30]}],
{n,2,20}]

**/