///////////////////////////////////////////////////////////////////////////////
version="$Id: ncalg.lib,v 1.11 2005-02-23 18:10:46 levandov Exp $";
category="Noncommutative";
info="
LIBRARY:  ncalg.lib      Definitions of important GR-algebras
AUTHORS:  Viktor Levandovskyy,     levandov@mathematik.uni-kl.de,
	  Oleksandr Motsak,        motsak@mathematik.uni-kl.de.

PROCEDURES:
sl(n[,p]);   returns U(sl_n) in char p, if an integer p is given,
sl2([p]);    returns U(sl_2) in the (e,f,h) presentation; in char p, if an integer p is given, 
g2([p]);     returns U(g_2) in the (x(i),y(i),Ha,Hb) presentation; in char p, if an integer p is given,
gl(n,[p]);    returns U(gl_n) in the (e_ij (1<i,j<n)) presentation; in char p, if an integer p is given,
Qso3([n]);   returns U_q(so_3) in the presentation of Klimyk, if integer n is given, the quantum parameter will be specialized at the (2n)-th root of unity, 
Qso3Casimir(n [,m]); returns a list with the (optionally normalized) Casimir elements of U_q(so_3) for the quantum parameter specialized at the n-th root of unity,
Qsl2([n]); returns ring, corresponding to the V_q(sl_2) and the ideal Qideal in it, such that U_q(sl_2) = V_q(sl_2) / Qideal; if n is specified, the quantum parameter q will be specialized at the n-th root of unity
Qsl3([n]); returns ring, corresponding to the V_q(sl_3) and the ideal Qideal in it, such that U_q(sl_3) = V_q(sl_3) / Qideal; if n is specified, the quantum parameter q will be specialized at the n-th root of unity
GKZsystem(A, sord, alg [,v]);  returns a ring and an ideal, describing
Gelfand-Kapranov-Zelevinsky system of differential equations
";

LIB "nctools.lib";
LIB "general.lib";
LIB "toric.lib"; // needed for GKZsystem

///////////////////////////////////////////////////////////////////////////////
static proc defInt ( list # )
// return 0 or int(#)
{
  int @p = 0;
  if ( size(#) > 0 )
  {
    if ( typeof( #[1] ) == "int" )
    {
      @p = #[1];
    };
  };
  return (@p);
};

///////////////////////////////////////////////////////////////////////////////
proc sl2(list #)
"USAGE:   sl2([p]), p an optional integer (field characteristic)
PURPOSE: set up the U(sl_2) in the (e,f,h) presentation over the field of char p
RETURN:  ring
NOTE:    you have to activate this ring with the 'setring' command
SEE ALSO: sl, g2, gl
EXAMPLE: example sl2; shows examples
"{
   int @p = defInt(#);
   ring @@@rrr=@p,(e,f,h),dp;
   matrix D[3][3]=0;
   D[1,2]=-h;
   D[1,3]=2*e;
   D[2,3]=-2*f;
   ncalgebra(1,D);
   return(@@@rrr);
}
example
{ "EXAMPLE:"; echo = 2;
   def a=sl2();
   setring a;
   a;
}

///////////////////////////////////////////////////////////////////////////////
proc sl(int n, list #)
"USAGE:   sl(n,[p]); n an integer, n>1; p an optional integer (field characteristic)
PURPOSE: set up the U(sl_n) in x(i),y(i),h(i) presentation over the field of char p
RETURN:  ring, describing 
NOTE:    You have to activate this ring with the 'setring' command. 
@*       The presentation of U(sl_n) is derived from the standard presentation of sl_n, 
@*       where positive resp. negative roots are denoted by x(i) resp. y(i) and
@*       Cartan elements are denoted by h(i).
SEE ALSO: sl2, g2, gl, Qsl3, Qso3
EXAMPLE: example sl; shows examples
"{
  if (n<2)
  {
    Print("Incorrect input");
    return(0);
  }
  if (n==2)
  {
    def @@@a=sl2(#);
    setring @@@a;
    return(@@@a);
  }
  
  int @p = defInt(#);
  ring @@@rr=@p,(x(1..n*(n-1)/2),y(1..n*(n-1)/2),h(1..n-1)),dp;
  intmat CNT[n][n]=0;
  matrix TMP[n][n]=0;
  int k,l=1,1;
  int buf=0;  
  list X,Y,H;
  for(k=1; k<=n; k++)
  {
    for(l=k+1; l<=n; l++)
    {
      buf = (l-k-1)*(2*n-l+k)/2 + k; 
      CNT[k,l] = buf;
      TMP[k,l] = 1;
      X[buf] = TMP;
      TMP = 0;
      CNT[l,k] = buf;
      TMP[l,k] = 1;
      Y[buf] = TMP;
      TMP=0;   
    }
  }
  for(k=1; k<=n-1; k++)
  {
    TMP[k,k]=1;
    TMP[k+1,k+1]=-1;
    H[k]=TMP;
    TMP=0;
  }  
  int i,j=1,1;
  number p,q=0,0;
  list V=X+Y+H;
  int v=size(V);
  matrix D[v][v]=0;  
  for(k=1; k<=v; k++)
  {
    for(l=k+1; l<=v; l++)
    {
      TMP=V[l]*V[k]-V[k]*V[l];
      for(i=1; i<=n; i++)
      {
        for(j=i+1; j<=n; j++)
        {
          buf=(j-i-1)*(2*n-j+i)/2+i;
          if (TMP[i,j]!=0)
          {
	    D[k,l]=D[k,l]+leadcoef(TMP[i,j])*x(buf);
          }
          if (TMP[j,i]!=0)
          {          
	    D[k,l]=D[k,l]+leadcoef(TMP[j,i])*y(buf);
          }
        }
      }
      i=1;
      while ( (TMP[i,i]==0) && (i<n) ) { i++; }
      for(j=i; j<=n-1; j++)
      {
        p=leadcoef(TMP[j,j]);
        q=leadcoef(TMP[j+1,j+1]);
        D[k,l]=D[k,l]+p*h(j);
	//        if ((j!=n-1)&&((p+q)!=0)) {D[k,l]=D[k,l]+(p+q)*h(j+1);}
        TMP[j+1,j+1]=TMP[j+1,j+1]+p; 
      }
    }
  }
  ncalgebra(1,D);
  return(@@@rr);
}
example
{ "EXAMPLE:"; echo = 2;
   def a=sl(3);
   setring a;
   a;
}

///////////////////////////////////////////////////////////////////////////////
proc g2(list #)
"USAGE:  g2([p]), p an optional integer (field characteristic)
PURPOSE: set up the U(g_2) in (x(i),y(i),Ha,Hb) presentation over the field of char p
RETURN:  ring
NOTE:    you have to activate this ring with the 'setring' command
SEE ALSO: sl, gl
EXAMPLE: example g2; shows examples
"{
  int @p = defInt(#);
  ring @@@rrr=@p,(x(1..6),y(1..6),Ha,Hb),dp;
  setring @@@rrr;
  matrix D[14][14];
  D[1,2]=-x(3);
  D[1,3]=-2*x(4);
  D[1,4]=3*x(5);
  D[1,7]=-Ha;
  D[1,9]=3*y(2);
  D[1,10]=2*y(3);
  D[1,11]=-y(4);
  D[1,13]=2*x(1);
  D[1,14]=-x(1);
  D[2,5]=x(6);
  D[2,8]=-Hb;
  D[2,9]=-y(1);
  D[2,12]=-y(5);
  D[2,13]=-3*x(2);
  D[2,14]=2*x(2);
  D[3,4]=3*x(6);
  D[3,7]=3*x(2);
  D[3,8]=-x(1);
  D[3,9]=-Ha-3*Hb;
  D[3,10]=-2*y(1);
  D[3,12]=-y(4);
  D[3,13]=-x(3);
  D[3,14]=x(3);
  D[4,7]=2*x(3);
  D[4,9]=-2*x(1);
  D[4,10]=-2*Ha-3*Hb;
  D[4,11]=y(1);
  D[4,12]=y(3);
  D[4,13]=x(4);
  D[5,7]=-x(4);
  D[5,10]=x(1);
  D[5,11]=-Ha-Hb;
  D[5,12]=y(2);
  D[5,13]=3*x(5);
  D[5,14]=-x(5);
  D[6,8]=-x(5);
  D[6,9]=-x(4);
  D[6,10]=x(3);
  D[6,11]=x(2);
  D[6,12]=-Ha-2*Hb;
  D[6,14]=x(6);
  D[7,8]=y(3);
  D[7,9]=2*y(4);
  D[7,10]=-3*y(5);
  D[7,13]=-2*y(1);
  D[7,14]=y(1);
  D[8,11]=-y(6);
  D[8,13]=3*y(2);
  D[8,14]=-2*y(2);
  D[9,10]=-3*y(6);
  D[9,13]=y(3);
  D[9,14]=-y(3);
  D[10,13]=-y(4);
  D[11,13]=-3*y(5);
  D[11,14]=y(5);
  D[12,14]=-y(6);
  ncalgebra(1,D);
  return(@@@rrr);
}
example
{ "EXAMPLE:"; echo = 2;
   def a = g2();
   setring a;  a;
}

///////////////////////////////////////////////////////////////////////////////
proc gl(int n, list #)
"
USAGE:   gl(n,[p]); n an integer, n>1;  p an optional integer (field characteristic)
PURPOSE: set up the U(gl_n) in the (e_ij (1<i,j<n)) presentation over the field of char p
RETURN:  ring
NOTE:    You have to activate this ring with the 'setring' command.
SEE ALSO: sl, g2
EXAMPLE: example gl; shows examples
"{
  if (n<2)
  {
    print("Incorrect input");
    return(0);
  };
  int @p = defInt(#);
  int i, j;
  string vs = "";
  for ( i = 1; i<= n ; i++ )
  {
        for ( j = 1; j<= n ; j++ )
	{
	    if ( vs != "" ) 
	    {
		vs = vs + ", ";
	    };
	    vs = vs + "e_" + string(i) + "_" + string(j);    
	};  
  };
  string strRING = "ring RING_GL=(" + string (@p) + "), (" + vs + "),dp;";
  execute( strRING );
  int N = nvars( RING_GL ); // n*n
  matrix D[N][N]=0;
  int k, l;
  int ik,il,jk,jl;
  poly p ;
  for( k=1; k<=N; k++)
  {
    ik = 1 + ((k-1)/n);
    jk = k -  n*(ik-1);
    
    for( l=k+1; l<=N; l++)
    {
        il = 1 + ((l-1)/n);
	jl = l -  n*(il-1);
	p = 0;    
	if( jl == ik )
	{
	    p = p + var ( (il-1)*n + jk );
	};
	if( jk == il )
	{
	    p = p - var ( (ik-1)*n + jl );
	};
	D[k,l]=p;
    };
  };
  ncalgebra(1,D);
  return(RING_GL);
}
example
{ "EXAMPLE:"; echo = 2;
   def a=gl(3);
   setring a; a;
};

///////////////////////////////////////////////////////////////////////////////
proc Qso3(list #)
"USAGE:   Qso3([n]), n an optional integers
PURPOSE: set up the U_q(so_3) in the presentation of Klimyk; if n is specified, the quantum parameter Q will be specialized at the (2*n)-th root of unity
RETURN:  ring
NOTE:    you have to activate this ring with the 'setring' command
SEE ALSO: sl, g2, gl, Qsl2, Qsl3, Qso3Casimir
EXAMPLE: example Qso3; shows examples
"{
  int @p = 2*defInt(#); 
  ring @@@r=(0,Q),(x,y,z),dp;
  minpoly = RootOfUnity(@p);
  matrix C[3][3];
  C[1,2]=Q2;
  C[1,3]=1/Q2;
  C[2,3]=Q2;
  matrix D[3][3];
  D[1,2]=-Q*z;
  D[1,3]=1/Q*y;
  D[2,3]=-Q*x;
  ncalgebra(C,D);
  return(@@@r);
}
example
{ "EXAMPLE:"; echo = 2;
   def K = Qso3(3);
   setring K;
   K;
}

///////////////////////////////////////////////////////////////////////////////
proc Qso3Casimir(int n, list #)
"USAGE:   Qso3Casimir(n [,m]), n an integer, m an optional integer
PURPOSE: compute the Casimir elements of U_q(so_3) for the quantum parameter specialized at the n-th root of unity; if m!=0 is given, polynomials will be normalized
RETURN:  list (of polynomials)
NOTE:    the result of the procedure makes sense only when the basering is U_q(so_3)
SEE ALSO: Qso3
EXAMPLE: example Qso3Casimir; shows examples
"{
  if ( npars(basering) !=1 )
  {
    "Error: wrong algebra. U_q(so3) has only one parameter";
    return(0);
  }
  if (n<1) { return(0); }
  number Q = par(1);
  int N=(n-1)/2;
  int NV=nvars(basering);
  number k1,k2;
  poly p,rs,hlp;
  list cp;
  int j;
  p=var(1);
  for(j=0; j<=N; j++)
  {
    k1 = binomial(n-j,j,0);
    k1=k1/(n-j);
    k1=k1*((-1)^j);
    k2=((Q^2)/(Q^4-1))^(2*j);
    k2=k2*k1;
    hlp=k2*(p)^(n-(2*j));
    rs=rs+hlp;
    hlp=0; k2=0; k1=0;
  }
  if (size(#)>0) 
  {
    int m = int(#[1]);
    if (m!=0)
    {
       rs = cleardenom(rs);
    }
  }
  cp[1] = rs;
  for(j=2; j<=NV; j++)
  {
    cp[j] = subst(rs,var(1),var(j));
  }
  return(cp);
}
example
{ "EXAMPLE:"; echo = 2;
   def R = Qso3(5);
   setring R;
   list C = Qso3Casimir(5);
   C;
   list Cnorm = Qso3Casimir(5,1);
   Cnorm;
}

///////////////////////////////////////////////////////////////////////////////
proc Qsl2(list #)
"USAGE:   Qsl2([n]), n an optional integer
PURPOSE: set up the U_q(sl_2) = V_q(sl_2) / Qideal via the ring V_q(sl_2) and the ideal 'Qideal' in it; 
@*       if n is specified, the quantum parameter q will be specialized at the n-th root of unity
RETURN:  ring (V_q(sl_2))
NOTE:    you have to activate this ring with the 'setring' command. 
@*       In order to create the U_q(sl_2) from the output, you have to call the command like 'qring Usl2q = Qideal;'
SEE ALSO: sl, Qsl3, Qso3 
EXAMPLE: example Qsl2; shows examples
"
{
  ring r=(0,q),(E,F,Ke,Kf),dp;
  int @p = defInt(#);
  if (@p >1) 
  {
    minpoly = RootOfUnity(@p);
  }
  matrix C = UpOneMatrix(4);;
  matrix D[4][4];
  C[1,3]=q^2;
  C[2,3]=1/(q^2);
  C[1,4]=1/(q^2);
  C[2,4]=q^2;
  D[1,2]=(1/(q-(1/q)))*(-Ke+Kf);
  ncalgebra(C,D);
  ideal Qideal = Ke*Kf-1;
  Qideal = twostd(Qideal);
  export Qideal; 
  return(r);
}
example
{ "EXAMPLE:"; echo = 2;
   def A = Qsl2(3);
   setring A;
   Qideal;
   qring Usl2q = Qideal;
   Usl2q;
}	

///////////////////////////////////////////////////////////////////////////////
proc Qsl3(list #)
"USAGE:   Qsl3([n]), n an optional integer
PURPOSE: set up the U_q(sl_3) = V_q(sl_3) / Qideal via the ring V_q(sl_3) and the ideal 'Qideal' in it; 
@*       if n is specified, the quantum parameter q will be specialized at the n-th root of unity
RETURN:  ring (V_q(sl_3))
NOTE:    you have to activate this ring with the 'setring' command. 
@*       In order to create the U_q(sl_3) from the output, you have to call the command like 'qring Usl3q = Qideal;'
SEE ALSO: sl, Qsl2, Qso3
EXAMPLE: example Qsl3; shows examples
"{
  int @p = defInt(#);
  ring @@@rrr=(0, q), (f12, f13, f23, k1, k2, l1, l2, e12, e13, e23), wp(2, 3, 2, 1, 1, 1, 1, 2, 3, 2);
  if (@p >1) 
  {
    minpoly = RootOfUnity(@p);
  }
  int @n = nvars(@@@rrr);
  matrix C = UpOneMatrix(@n);
  matrix D[@n][@n];
  // some constants
  number q1 =    1/q;
  number Q  = (q )^2;
  number Q1 = (q1)^2;
  //   number QQ = Q - Q1; // q2 - 1/(q2)
  number QQ1= 1 / (Q - Q1);  
  // relations:
  C[1,2] = Q1;
  C[2,3] = C[1,2];
  C[8,9] = C[1,2];
  C[9,10]= C[1,2];
  C[1,3] = Q;
  C[8,10]= C[1,3];
  
  D[1,3] = -q*(f13);
  D[8,10]= -q*(e13);
  // V_q(sl_3)
  D[1,8] = QQ1 * ( (k1) ^ 2 - (l1) ^ 2 );
  D[3,10]= QQ1 * ( (k2) ^ 2 - (l2) ^ 2 );
  D[2,9] = -QQ1 * ( ((k1)^2)*((k2)^2) - ((l1)^2)*((l2)^2) );
  D[2, 8] =   q * (f23) * ((k1)^2);
  D[3, 9] =   q * ((k2)^2) * (e12);
  D[1, 9] = -q1 * ((l1)^2) * (e23);
  D[2, 10]= -q1 * (f12) * ((l2)^2);
  // k1
  C[ 4, 8 ]= Q1;
  C[ 4, 9 ]= q1;
  C[ 4, 10]= q;
  // l1
  C[ 6, 8 ]= Q;
  C[ 6, 9 ]= q;
  C[ 6, 10]= q1;
  // k2
  C[ 5, 8 ]= q;
  C[ 5, 9 ]= q1;
  C[ 5, 10]= Q1;
  // l2
  C[ 7, 8 ]= q1;
  C[ 7, 9 ]= q;
  C[ 7, 10]= Q;
  // k1
  C[ 1, 4 ]= Q1;
  C[ 2, 4 ]= q1;
  C[ 3, 4 ]= q;
  // l1
  C[ 1, 6 ]= Q;
  C[ 2, 6 ]= q;
  C[ 3, 6 ]= q1;
  // k2
  C[ 1, 5 ]= q;
  C[ 2, 5 ]= q1;
  C[ 3, 5 ]= Q1;
  // l2
  C[ 1, 7 ]= q1;
  C[ 2, 7 ]= q;
  C[ 3, 7 ]= Q;
  ncalgebra(C,D); // the V_q(sl3) is done
  ideal Qideal = k1*l1-1,  k2*l2-1;
  Qideal = twostd(Qideal);
  export Qideal;
  return(@@@rrr);
}
example
{ "EXAMPLE:"; echo = 2;
   def B = Qsl3(5);
   setring B;
   qring Usl3q = Qideal;
   Usl3q;
}

proc GKZsystem(intmat A, string sord, string alg, list #)
"USAGE:   GKZsystem(A, sord, alg, [,v]), where 
@*        A    is an intmat, defining the system, 
@*        sord is a string with desired term ordering, 
@*        alg  is a string, saying which algorithm to use (exactly like in toric_lib), 
@*        v    is an optional intvec.
PURPOSE: compute the ring and the GKZ system of equations in it
RETURN:  ring (moreover, the ideal GKZid with equations is exported in it)
NOTE:    you have to activate the ring with the 'setring' command; 
the procedure is elaborated by Oleksandr Yena
OVERVIEW: this procedure uses toric_lib and therefore inherits its input requirements: 
@*        possible values for input variable 'alg' are: \"ect\",\"pt\",\"blr\", \"hs\", \"du\". 
@*        As for the term ordering, it should be a string in Singular format like \"lp\",\"dp\", etc. @*        Please consult with the toric_lib for allowed orderings and more details.
SEE ALSO: toric_lib
EXAMPLE: example GKZsystem; shows examples
"{
  int d = nrows(A);
  int n = ncols(A);
  execute("ring r1=0,(d(1..n)),"+sord+";");
  ideal I0;
  if (size(#)==0)
  {
    I0 = toric_ideal(A, alg);
  }
  else 
  {
    if ( typeof(#[1]) == "intvec" )
    {
      intvec V = intvec(#[1]);
      I0 = toric_ideal(A, alg, V);
    }
    else
    {
      "Wrong type of the optional argument. Intvec expected.";
      return();
    }
  };
  execute("ring GR = (0,b(1..d)),(x(1..n),d(1..n)),"+sord+";");
  Weyl();
  int i,j;
  poly p;
  ideal I;
  for (i=1; i<=d; i++)
  {
    p = -b(i);
    for (j=1; j<=n; j++)
    {
      p = p+ A[i,j]*x(j)*d(j);
    }
  I = I, p;
  }
  I = I, imap(r1,I0);
  I = simplify(I,2);
  ideal GKZid = I;
  export GKZid;
  return(GR);
}
example
{"EXAMPLE:"; echo = 2;
  intmat A[2][4]=3,2,1,0,0,1,2,3;
  print(A);
  def D1 = GKZsystem(A,"lp","ect");
  setring D1;
  D1;
  print(GKZid);
  intvec v=1,1,1,1;
  def D2 = GKZsystem(A,"lp","blr",v);
  setring D2;
  print(GKZid);
}

///////////////////////////////////////////////////////////////////////////////