///////////////////////////////////////////////////////////////////////////////
version="$Id: ncdecomp.lib,v 1.13 2007-11-28 16:36:46 motsak Exp $";
category="Noncommutative";
info="
LIBRARY:  ncdecomp.lib     Decomposition of a module into its central characters
AUTHORS:  Viktor Levandovskyy,     levandov@mathematik.uni-kl.de.

OVERVIEW:
This library presents algorithms for the  central character
decomposition of a module, i.e. a decomposition into generalized weight modules with respect to the center. Based on ideas of O. Khomenko and V. Levandovskyy (see the article [L2] in the References for details).

PROCEDURES:
CentralQuot(M,G);       central quotient M:G,
CentralSaturation(M,T); central saturation ((M:T):...):T) ( = M:T^{\infty}),
CenCharDec(I,C);        decomposition of I into central characters w.r.t. C
IntersectWithSub(M,Z);  intersection of M with the subalgebra, generated by pairwise commutative elements of Z.
";

  LIB "ncalg.lib";
  LIB "primdec.lib";
///////////////////////////////////////////////////////////////////////////////
static proc CharKernel(list L, int i)
{
// compute \cup L[j], j!=i
  int sL = size(L);
  if ( (i<=0) || (i>sL))  { return(0); }
  int j;
  list Li;
  if (i ==1 )
  {
    Li = L[2..sL];
  }
  if (i ==sL )
  {
    Li = L[1..sL-1];
  }
  if ( (i>1) && (i < sL))
  {
    Li = L[1..i-1];
    for (j=i+1; j<=sL; j++)
    {
      Li[j-1] = L[j];
    }
  }
//  print("intersecting kernels...");
  module Cres = intersect(Li[1..size(Li)]);
  return(Cres);
}
///////////////////////////////////////////////////////////////////////////////
static proc CentralQuotPoly(module M, poly g)
{
// here an elimination of components should be used !
  int N=nrows(M); // M = A^N /I_M
  module @M;
  int i,j;
  for(i=1; i<=N; i++)
  {
   @M=@M,g*gen(i);
  }
  @M = simplify(@M,2);
  @M = @M,M;
  module S = syz(@M);
  matrix s = S;
  module T;
  vector t;
  for(i=1; i<=ncols(s); i++)
  {
    t = 0*gen(N);
    for(j=1; j<=N; j++)
    {
      t = t + s[j,i]*gen(j);
    }
    T[i] = t;
  }
  T = simplify(T,2);
  return(T);
}
///////////////////////////////////////////////////////////////////////////////
static proc MyIsEqual(module A, module B)
{
// both A and B are submodules of free module
  option(redSB);
  option(redTail);
  if (attrib(A,"isSB")!=1)
  {
    A = std(A);
  }
  if (attrib(B,"isSB")!=1)
  {
    B = std(B);
  }
  int ANSWER = 1;
  if ( ( ncols(A) == ncols(B) ) && ( nrows(A) == nrows(B) ) )
  {
    module @AB = A-B;
    @AB = simplify(@AB,2);
    if (@AB[1]!=0) { ANSWER = 0; }
  }
  else { ANSWER = 0; }
  return(ANSWER);
}
///////////////////////////////////////////////////////////////////////////////
proc CentralQuot(module I, ideal G)
"USAGE:  CentralQuot(M, G), M a module, G an ideal
ASSUME: G is an ideal in the center
RETURN:  module
PURPOSE: compute the central quotient M:G
THEORY:  for an ideal G of the center of an algebra and a submodule M of A^n, the central quotient of M by G is defined to be
@* M:G  :=  { v in A^n | z*v in M, for all z in G }.
NOTE:    the output module is not necessarily given in a Groebner basis
SEE ALSO: CentralSaturation, CenCharDec
EXAMPLE: example CentralQuot; shows examples
"{
  int i;
  list @L;
  for(i=1; i<=size(G); i++)
  {
    @L[i] = CentralQuotPoly(I,G[i]);
  }
  module @I = intersect(@L[1..size(G)]);
  if (nrows(@I)==1)
  {
    @I = ideal(@I);
  }
  return(@I);
}
example
{ "EXAMPLE:"; echo = 2;
   option(returnSB);
   def a = makeUsl2();
   setring a;
   ideal I = e3,f3,h3-4*h;
   I = std(I);
   poly C=4*e*f+h^2-2*h;
   ideal G = (C-8)*(C-24);
   ideal R = CentralQuot(I,G);
   R;
}
///////////////////////////////////////////////////////////////////////////////
proc CentralSaturation(module M, ideal T)
"USAGE:  CentralSaturation(M, T), for a module M and an ideal T
ASSUME: T is an ideal in the center
RETURN:  module
PURPOSE: compute the central saturation of M by T, that is M:T^{\infty}, by repititive application of @code{CentralQuot}
NOTE:    the output module is not necessarily a Groebner basis
SEE ALSO: CentralQuot, CenCharDec
EXAMPLE: example CentralSaturation; shows examples
"{
  option(redSB);
  option(redTail);
  option(returnSB);
  module Q=0;
  module S=M;
  while ( !MyIsEqual(Q,S) )
  {
    Q = CentralQuot(S, T);
    S = CentralQuot(Q, T);
  }
  if (nrows(Q)==1)
  {
    Q = ideal(Q);
  }
//  Q = std(Q);
  return(Q);
}
example
{ "EXAMPLE:"; echo = 2;
   option(returnSB);
   def a = makeUsl2();
   setring a;
   ideal I = e3,f3,h3-4*h;
   I = std(I);
   poly C=4*e*f+h^2-2*h;
   ideal G = C*(C-8);
   ideal R = CentralSaturation(I,G);
   R=std(R);
   vdim(R);
   R;
}
///////////////////////////////////////////////////////////////////////////////
proc CenCharDec(module I, def #)
"USAGE:  CenCharDec(I, C);  I a module, C an ideal
ASSUME: C consists of generators of the center
RETURN:  a list L, where each entry consists of three records (if a finite decomposition exists)
@*       L[*][1] ('ideal' type), the central character as the maximal ideal in the center,
@*       L[*][2] ('module' type), the Groebner basis of the weight module, corresponding to the character in  L[*][1],
@*       L[*][3] ('int' type) is the vector space dimension of the weight module (-1 in case of infinite dimension);
PURPOSE: compute a finite decomposition of C into central characters or determine that there is no finite decomposition
NOTE:     actual decomposition is a sum of L[i][2] above;
@*        some modules have no finite decomposition (in such case one gets warning message)
SEE ALSO: CentralQuot, CentralSaturation
EXAMPLE: example CenCharDec; shows examples
"
{
  list Center;
  if (typeof(#) == "ideal")
  {
    int cc;
    ideal tmp = ideal(#);
    for (cc=1; cc<=size(tmp); cc++)
    {
      Center[cc] = tmp[cc];
    }
    kill tmp;
  }
  if (typeof(#) == "list")
  {
    Center = #;
  }
// M = A/I
//1. Find the Zariski closure of Supp_Z M
// J = Ann_M 1 == I
// J \cap Z:
  option(redSB);
  option(redTail);
  option(returnSB);
  def @A = basering;
  setring @A;
  int sZ=size(Center);
  int i,j;
  poly t=1;
  for(i=1; i<=nvars(@A); i++)
  {
    t=t*var(i);
  }
  ring @Z=0,(@z(1..sZ)),dp;
//  @Z;
  def @ZplusA = @A+@Z;
  setring @ZplusA;
//  @ZplusA;
  ideal I     = imap(@A,I);
  list Center = imap(@A,Center);
  poly t      = imap(@A,t);
  ideal @Ker;
  for(i=1; i<=sZ; i++)
  {
    @Ker[i]=@z(i) - Center[i];
  }
  @Ker = @Ker,I;
  ideal @JcapZ = eliminate(@Ker,t);
// do not forget parameters of a basering!
  string strZ="ring @@Z=("+charstr(@A)+"),(@z(1.."+string(sZ)+")),dp;";
//  print(strZ);
  execute(strZ);
  setring @@Z;
  ideal @JcapZ = imap(@ZplusA,@JcapZ);
  @JcapZ = std(@JcapZ);
//  @JcapZ;
  int sJ = vdim(@JcapZ);
  if (sJ==-1)
  {
    "There is no finite decomposition";
    return(0);
  }
//  print(@JcapZ);
// 2. compute the min.ass.primes of the ideal in the center
  list @L = minAssGTZ(@JcapZ);
  int sL = size(@L);
//  print("etL:");
//  @L;
// exception: is sL==1, the whole ideal has unique cen.char
  if (sL ==1)
  {
    setring @A;
    map @M = @@Z,Center[1..size(Center)];
    list L = @M(@L);
    list @R;
    @R[1] = L[1];
    if (nrows(@R[1])==1)
    {
      @R[1] = ideal(@R[1]);
    }
    @R[2] = I;
    if (nrows(@R[2])==1)
    {
      @R[2] = ideal(@R[2]);
    }
    @R[2] = std(@R[2]);
    @R[3] = vdim(@R[2]);
    return(@R);
  }
  list @CharKer;
  for(i=1; i<=sL; i++)
  {
    @L[i] = std(@L[i]);
  }
// 3. compute the intersections of characters
  for(i=1; i<=sL; i++)
  {
    @CharKer[i] = CharKernel(@L,i);
  }
//  print("Charker:");
//  @CharKer;
// 4. Go back to the algebra and compute central saturations
  setring @A;
  map @M = @@Z,Center[1..size(Center)];
  list L = @M(@CharKer);
  list R,@R;
  for(i=1; i<=sL; i++)
  {
    @R[1] = L[i];
    if (nrows(@R[1])==1)
    {
      @R[1] = ideal(@R[1]);
    }
    @R[2] = CentralSaturation(I,L[i]);
    if (nrows(@R[2])==1)
    {
      @R[2] = ideal(@R[2]);
    }
    @R[2] = std(@R[2]);
    @R[3] = vdim(@R[2]);
     R[i] = @R;
  }
  return(R);
}
example
{ "EXAMPLE:"; echo = 2;
   option(returnSB);
   def a = makeUsl2(); // U(sl_2) in characteristic 0
   setring a;
   ideal I = e3,f3,h3-4*h;
   I = twostd(I);           // two-sided ideal generated by I
   vdim(I);                 // it is finite-dimensional
   ideal Cn = 4*e*f+h^2-2*h; // the only central element
   list T = CenCharDec(I,Cn);
   T;
   // consider another example
   ideal J = e*f*h;
   CenCharDec(J,Cn);
}
///////////////////////////////////////////////////////////////////////////////
proc IntersectWithSub (ideal M, def #)
"USAGE:  IntersectWithSub(M,Z),  M an ideal, Z an ideal
ASSUME: Z consists of pairwise commutative elements
RETURN:  ideal, of two-sided generators, not a Groebner basis
PURPOSE: computes an intersection of M with the subalgebra, generated by Z
NOTE:    usually Z consists of generators of the center
EXAMPLE: example IntersectWithSub; shows an example
"
{
  ideal Z;
  if (typeof(#) == "list")
  {
    int cc;
    list tmp = #;
    for (cc=1; cc<=size(tmp); cc++)
    {
      Z[cc] = tmp[cc];
    }
    kill tmp;
  }
  if (typeof(#) == "ideal")
  {
    Z = #;
  }
  // returns a submodule of M, equal to M \cap Z
  // correctness: Z should consists of pairwise
  // commutative elements
  int nz = size(Z);
  int i,j;
  poly p;
  for (i=1; i<nz; i++)
  {
    for (j=i+1; j<=nz; j++)
    {
      p = bracket(Z[i],Z[j]);
      if (p!=0)
      {
        "Error: generators of the subalgebra do not commute.";
        return(ideal(0));
      }
    }
  }
  // main action
  def B = basering;
  setring B;
  string s1,s2;
  s1 = "ring @Z = (";
  s2 = s1 + charstr(basering) + "),(z(1.." + string(nz)+")),Dp";
  //  s2;
  execute(s2);
  setring B;
  map F = @Z,Z;
  setring @Z;
  ideal PreM = preimage(B,F,M);
  PreM = std(PreM);
  setring B;
  ideal T = F(PreM);
  return(T);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring R=(0,a),(e,f,h),Dp;
  matrix @d[3][3];
  @d[1,2]=-h;
  @d[1,3]=2e;
  @d[2,3]=-2f;
  def r = nc_algebra(1,@d); setring r; // parametric U(sl_2)
  ideal I = e,h-a;
  ideal C;
  C[1] = h^2-2*h+4*e*f; // the center of U(sl_2)
  ideal X = IntersectWithSub(I,C);
  X;
  ideal G = e*f, h; // the biggest comm. subalgebra of U(sl_2)
  ideal Y = IntersectWithSub(I,G);
  Y;
}