///////////////////////////////////////////////////////////////////////////////
version="$Id$";
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^infinity),
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";
  LIB "central.lib";

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

proc testncdecomplib()
{
  example CentralQuot;
  example CentralSaturation;
  example CenCharDec;
  example IntersectWithSub;
}

static proc CharKernel(list L, int i)
{
 // todo: think on more effective way of doing it...
// 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)]); // uses std, try modulo!
  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 = slimgb(A);
  }
  if (attrib(B,"isSB")!=1)
  {
    B = slimgb(B);
  }
  int ANSWER = 1;
  if ( ( ncols(A) == ncols(B) ) && ( nrows(A) == nrows(B) ) )
  {
    module @AB = module(matrix(A)-matrix(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 of the base ring
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
"{
/* check assupmtion. Elt's of G must be central */
  if (! inCenter(G) )
  {
    ERROR("ideal in the 2nd argument is not in the center of the base ring!");
  }
  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;  // C in Z(U(sl2)), the central element
   ideal G = (C-8)*(C-24);  // G normal factor in Z(U(sl2)) as an ideal in the center
   ideal R = CentralQuot(I,G);  // same as 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 of the base ring
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
"{
/* check assupmtion. Elt's of T must be central */
  if (! inCenter(T) )
  {
    ERROR("ideal in the 2nd argument is not in the center of the base ring!");
  }
  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 of the base ring
RETURN:  a list L, where each entry consists of three records (if a finite decomposition exists)
@*       L[*][1] ('ideal' type), the central character as a 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 the sum of L[i][2] above;
@*        some modules have no finite decomposition (in such case one gets warning message)
@*        The function @code{central} in @code{central.lib} may be used to obtain C, when needed.
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 = #;
  }

/* check assupmtion. Elt's of G must be central */
  if (! inCenter(Center) )
  {
    ERROR("ideal in the 2nd argument is not in the center of the base ring!");
  }
  int ppl = printlevel-voice+2;
// 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);
  dbprint(ppl,"// -1-1- starting the computation of preimage in Z");
  dbprint(ppl-1, @Ker);
  ideal @JcapZ = slimgb(@Ker);
  @JcapZ = nselect(@JcapZ,intvec(1..nvars(@A)));
  dbprint(ppl,"// -1-2- finished the computation of preimage in Z");
  dbprint(ppl-1, @JcapZ);
// do not forget parameters of a basering!
// hmmm: todo ringlist
  string strZ="ring @@Z=("+charstr(@A)+"),(@z(1.."+string(sZ)+")),dp;";
//  print(strZ);
  execute(strZ);
  setring @@Z;
  ideal @JcapZ = imap(@ZplusA,@JcapZ);
  dbprint(ppl,"// -1-3- starting the cosmetic Groebner basis in Z");
  @JcapZ = slimgb(@JcapZ); // evtl. groebner?
//  @JcapZ;
  dbprint(ppl,"// -1-4- finished the cosmetic Groebner basis in Z");
  dbprint(ppl-1, @JcapZ);
  int sJ = vdim(@JcapZ);
  dbprint(ppl,"// -1-5- the K-dimension of support is "+string(sJ));
  if (sJ==-1)
  {
    "There is no finite decomposition";
    return(0);
  }
//  print(@JcapZ);
// 2. compute the min.ass.primes of the ideal in the center
  dbprint(ppl,"// -2-1- starting the computation of minimal primes in Z");
  list @L = minAssGTZ(@JcapZ);
  int sL = size(@L);
  dbprint(ppl,"// -2-2- finished the computation of " + string(sL)+ " minimal primes in Z");
//  print("etL:");
//  @L;
// exception: is sL==1, the whole ideal has unique cen.char
  if (sL ==1)
  {
    dbprint(ppl-1,"// -2-3- the whole module is gen. weight module itself");
    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]);
    }
    dbprint(ppl-1,"// -2-4- final cosmetic Groebner basis");
    @R[2] = slimgb(@R[2]);
    @R[3] = vdim(@R[2]);
    return(list(@R)); // for compliance with output a list
  }
  dbprint(ppl-1,"// -2-3- there are several characters");
  dbprint(ppl,"// -*- computing Groebner bases of components (commutative)");
  list @CharKer;
  for(i=1; i<=sL; i++)
  {
    @L[i] = slimgb(@L[i]);
  }
  dbprint(ppl,"// -*- finished computing Groebner bases of components");
// 3. compute the intersections of characters
  dbprint(ppl,"// -3- compute the intersections of characters");
  for(i=1; i<=sL; i++)
  {
    @CharKer[i] = CharKernel(@L,i);
  }
  dbprint(ppl,"// -3- the intersections of characters is done");
  //  dbprint(ppl-1,@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;
  dbprint(ppl,"// -4- compute the central saturations");
  dbprint(ppl-1,L);
  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] = slimgb(@R[2]);
    @R[3] = vdim(@R[2]);
     R[i] = @R;
  }
  dbprint(ppl,"// -4- central saturations are done");
  return(R);
}
example
{ "EXAMPLE:"; echo = 2; printlevel=0;
   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 the intersection of M with the subalgebra, generated by Z
NOTE:    usually Z consists of generators of the center
@* The function @code{central} from @code{central.lib} may be used to obtain the center Z, if needed.
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
  // assume/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;
  // todo: make ringlist from it!
  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); // reformulate using gb engine? todo?
  PreM = slimgb(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;
}