source: git/Singular/LIB/ratgb.lib @ 380a17b

/////////////////////////////////////////////////////////////////////////
version="version ratgb.lib 4.0.0.0 Jun_2013 ";
category="Noncommutative";
info="
LIBRARY: ratgb.lib  Groebner bases in Ore localizations of noncommutative G-algebras
AUTHOR: Viktor Levandovskyy,     levandov@risc.uni-linz.ac.at

OVERVIEW:
Theory: Let A be an operator algebra with @code{R = K[x1,.,xN]} as subring.
The operators are usually denoted by @code{d1,..,dM}.

Assume, that A is a @code{G}-algebra, then the set @code{S=R-0} is multiplicatively
closed Ore set in A.
That is, for any s in S and a in A, there exist t in S and b in A, such that @code{sa=bt}.
In other words, one can transform any left fraction into a right fraction.
The algebra @code{A_S} is called an Ore localization of A with respect to S.

This library provides Groebner basis procedure for A_S, performing polynomial (that is
fraction-free) computations only. Note, that there is ongoing development of the
subsystem called Singular:Locapal, which will provide yet another approach to Groebner
bases over such Ore localizations.

Assumptions: in order to treat such localizations constructively, some care need to be taken.
We will assume that the variables @code{x1,...,xN} from above (which will become invertible
in the localization) come as the first block among the variables of the basering.
Moreover, the ordering on the basering must be an antiblock ordering, that is its
matrix form has the left upper @code{NxN} block zero. Here is a recipe to create such
an ordering easily: use 'a(w)' definitions of the ordering N times with intvecs @code{w_i}
of the following form: @code{w_i} has first N components zero. The rest entries need to
be positive and such, that @code{w1,..,wN} are linearly independent (see an example below).

Guide: with this library, it is possible
- to compute a Groebner basis of an ideal or a submodule in the 'rational'
  Ore localization D = A_S
- to compute a dimension of associated graded submodule (called D-dimension)
- to compute a vector space dimension over Quot(R) of a submodule of
  D-dimension 0 (so called D-finite submodule)
- to compute a basis over Quot(R) of a D-finite submodule

PROCEDURES:
ratstd(I, n); compute Groebner basis and dimensions in Ore localization of the basering

Support: SpezialForschungsBereich F1301 of the Austrian FWF and
Transnational Access Program of RISC Linz, Austria

SEE ALSO: jacobson_lib
";

LIB "poly.lib";
LIB "dmodapp.lib"; // for Appel1, Appel2, Appel4


static proc rm_content_id(def J)
"USAGE:  rm_content_id(I);  I an ideal/module
RETURN:  the same type as input
PURPOSE: remove the content of every generator of I
EXAMPLE: example rm_content_id; shows examples
"
{
  def  I = J;
  int i;
  int s = ncols(I);
  for (i=1; i<=s; i++)
  {
    if (I[i]!=0)
    {
      I[i] = I[i]/content(I[i]);
    }
  }
  return(I);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = (0,k,n),(K,N),dp;
  ideal I = n*((k+1)*K - (n-k)), k*((n-k+1)*N - (n+1));
  I;
  rm_content_id(I);
  module M = I[1]*gen(1), I[2]*gen(2);
  print(rm_content_id(M));
}

proc ratstd(def I, int is, list #)
"USAGE:  ratstd(I, n [,eng]);  I an ideal/module, n an integer, eng an optional integer
RETURN:  ring
PURPOSE: compute the Groebner basis of I in the Ore localization of
the basering with respect to the subalgebra, generated by first n variables
ASSUME: the variables of basering are organized in two blocks and
- the first block of length n contains the elements with respect to which one localizes,
- the basering is equipped with the elimination block ordering for the variables
  in the second block
NOTE: the output ring C is commutative. The ideal @code{rGBid} in C
represents the rational form of the output ideal @code{pGBid} in the basering.
- During the computation, the D-dimension of I and the corresponding
dimension as K(x)-vector space of I are computed and printed out.
- Setting optional integer eng to 1, slimgb is taken as Groebner engine
DISPLAY: In order to see the steps of the computation, set printlevel to >=2
EXAMPLE: example ratstd; shows examples
"
{
  int ppl = printlevel-voice+1;
  int eng = 0;
  // optional arguments
  if (size(#)>0)
  {
    if (typeof(#[1]) == "int")
    {
      eng = int(#[1]);
    }
  }

  dbprint(ppl,"engine chosen to be");
  dbprint(ppl,eng);

  // 0. do the subst's /reformulations
  // for the time being, ASSUME
  // the ord. is an elim. ord. for D
  // and the block of X's is on the left
  // its length is 'is'

  int i,j,k;
  dbprint(ppl,"// -1- creating K(x)[D]");

  // 1. create K(x)[D], commutative
  def save = basering;
  list L = ringlist(save);
  list RL, tmp1,tmp2,tmp3,tmp4;
  intvec iv;
  // copy: field, enlarge it with Xs

  if ( size(L[1]) == 0)
  {
    // i.e. the field with char only
    tmp2[1] = L[1];
    //    tmp1 = L[2];
    j    = size(L[2]);
    iv   = 1;
    for (i=1; i<=is; i++)
    {
      tmp1[i] = L[2][i];
      iv = iv,1;
    }
    iv = iv[1..size(iv)-1]; //extra 1
    tmp2[2]       = tmp1;
    tmp3[1] = "lp";
    tmp3[2] = iv;
    //    tmp2[3] = 0;
    tmp4[1] = tmp3;
    tmp2[3] = tmp4;
    //[1] = "lp";
    //    tmp2[3][2] = iv;
    tmp2[4]       = ideal(0);
    RL[1] = tmp2;
  }

  if ( size(L[1]) >0 )
  {
    // TODO!!!!!
    tmp2[1] = L[1][1]; //char K
    // there are parameters
    // add to them X's, IGNORE alg.extension
    // the ordering on pars
    tmp1 = L[1][2]; // param names
    j    = size(tmp1);
    iv = L[1][3][1][2];
    for (i=1; i<=is; i++)
    {
      tmp1[j+i] = L[2][i];
      iv = iv,1;
    }
    tmp2[2] = tmp1;
    tmp2[3] = L[1][3];
    tmp2[3][1][2] = iv;
    tmp2[4] = ideal(0);
    RL[1] = tmp2;
  }

  // vars: leave only D's
  kill tmp1; list tmp1;
  //  tmp1 = L[2];
  for (i=is+1; i<= size(L[2]); i++)
  {
    tmp1[i-is] = L[2][i];
  }
  RL[2] = tmp1;

  // old: assume the ordering is the block with (a(0:is),ORD)
  // old :set up ORD as the ordering
  //  L; "RL:"; RL;
  if (size(L[3]) != 3)
  {
    //"note: strange ordering";
    // NEW assume: ordering is the antiblock with (a(0:is),a(*1),a(*), ORD)
    // get the a() parts after is => they should form a complete D-ordering
    list L3 = L[3]; list NL3; kill tmp3; list tmp3;
    int @sl = size(L3);
    int w=1; int z; intvec va,vb;
    while(L3[w][1] == "a")
    {
      va = L3[w][2];
      for(z=1;z<=nvars(save)-is;z++)
      {
        vb[z] = va[is+z];
      }
      tmp3[1] = "a";
      tmp3[2] = vb;
      NL3[w] = tmp3; tmp3=0;
      w++;
    }
    // check for completeness: must be >= nvars(save)-is rows
    if (w < nvars(save)-is)
    {
      "note: ordering is incomplete on D. Adding lower Dp block";
      // adding: positive things like Dp
      tmp3[1]= "Dp";
      for (z=1; z<=nvars(save)-is; z++)
      {
        va[is+z] = 1;
      }
      tmp3[2] = va;
      NL3[w] = tmp3; tmp3 = 0;
      w++;
    }
    NL3[w] = L3[@sl]; // module ord?
    RL[3] = NL3;
  }
  else
  {
    kill tmp2; list tmp2;
    tmp2[1] = L[3][2];
    tmp2[2] = L[3][3];
    RL[3]   = tmp2;
  }
  // factor ideal is ignored
  RL[4] = ideal(0);

  //  "ringlist:"; RL;

  def @RAT = ring(RL);

  dbprint(ppl,"// -2- preprocessing with content");
  // 2. preprocess input with rm_content_id
  setring @RAT;
  dbprint(ppl-1, @RAT);
  //  ideal CI = imap(save,I);
  def CI = imap(save,I);
  CI = rm_content_id(CI);
  dbprint(ppl-1, CI);

  dbprint(ppl,"// -3- running groebner");
  // 3. compute G = GB(I) w.r.t. the elim. ord. for D
  setring save;
  //  ideal CI = imap(@RAT,CI);
  def CI = imap(@RAT,CI);
  option(redSB);
  option(redTail);
  if (eng)
  {
    def G = slimgb(CI);
  }
  else
  {
    def G = groebner(CI);
  }
  //  ideal G = groebner(CI); // although slimgb looks better
  // def G = slimgb(CI);
  G = simplify(G,2); // to be sure there are no 0's
  dbprint(ppl-1, G);

  dbprint(ppl,"// -4- postprocessing with content");
  // 4. postprocess the output with 1) rm_content_id,  2) lm-minimization;
  setring @RAT;
  // ideal CG = imap(save,G);
  def CG = imap(save,G);
  CG = rm_content_id(CG);
  CG = simplify(CG,2);
  dbprint(ppl-1, CG);

  // warning: a bugfarm! in this ring, the ordering might change!!! (see appelF4)
  // so, simplify(32) should take place in the orig. ring! and NOT here
  //  CG = simplify(CG,2+32);

  // 4b. create L(G) with X's as coeffs (for minimization)
  setring save;
  G = imap(@RAT,CG);
  int sG  = ncols(G);
  //  ideal LG;
  def LG = G;
   for (i=1; i<= sG; i++)
   {
     LG[i] = lead(G[i]);
   }
  // compute the D-dimension of the ideal in the ring @RAT
  setring @RAT;
  //  ideal LG = imap(save,LG);
  def LG = imap(save,LG);
  //  ideal LGG = groebner(LG); // cosmetics
  def LGG = groebner(LG); // cosmetics
  int d = dim(LGG);
  int Ddim = d;
  printf("the D-dimension is %s",d);
  if (d==0)
  {
    d = vdim(LGG);
    int Dvdim = d;
    printf("the K-dimension is %s",d);
  }
  //  ideal SLG = simplify(LG,8+32); //contains zeros
  def SLG = simplify(LG,8+32); //contains zeros
  setring save;
  //  ideal SLG = imap(@RAT,SLG);
  def SLG = imap(@RAT,SLG);
  // simplify(LG,8+32); //contains zeros
  intvec islg;
  if (SLG[1] == 0)
  {  islg = 0;  }
  else
  {    islg = 1;  }
  for (i=2; i<= ncols(SLG); i++)
  {
    if (SLG[i] == 0)
    {
      islg = islg, 0;
    }
    else
    {
      islg = islg, 1;
    }
  }
  for (i=1; i<= ncols(LG); i++)
  {
    if (islg[i] == 0)
    {
      G[i] = 0;
    }
  }
  G = simplify(G,2); // ready!
  //  G = imap(@RAT,CG);
  // return the result
  //  ideal pGBid = G;
  def pGBid = G;
  export pGBid;
  //  export Ddim;
  //  export Dvdim;
  setring @RAT;
  //  ideal rGBid = imap(save,G);
  def rGBid = imap(save,G);
  // CG;
  export rGBid;
  setring save;
  return(@RAT);
  //  kill @RAT;
  //  return(G);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = (0,c),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
  // this ordering is an antiblock ordering, as it must be
  def S = Weyl(); setring S;
  // the ideal I below annihilates parametric Appel F4 function
  // where we set parameters to a=-2, b=-1 and d=0
  ideal I =
    x*Dx*(x*Dx+c-1) - x*(x*Dx+y*Dy-2)*(x*Dx+y*Dy-1),
    y*Dy*(y*Dy-1) - y*(x*Dx+y*Dy-2)*(x*Dx+y*Dy-1);
  int is = 2; // hence 1st and 2nd variables, that is x and y
  // will become invertible in the localization
  def A = ratstd(I,2); // main call
  pGBid; // polynomial form of the basis in the localized ring
  setring A; // A is a commutative ring used for presentation
  rGBid; // "rational" or "localized" form of the basis
  //--- Now, let us compute a K(x,y) basis explicitly
  print(matrix(kbase(rGBid)));
}

/*
oldExampleForDoc()
{
  // VL: removed since it's too easy
  ring r = 0,(k,n,K,N),(a(0,0,1,1),a(0,0,1,0),dp);
  // note, that the ordering must be an antiblock ordering
  matrix D[4][4]; D[1,3] = K; D[2,4] = N;
  def S = nc_algebra(1,D);
  setring S; // S is the 2nd shift algebra
  ideal I = (k+1)*K - (n-k), (n-k+1)*N - (n+1);
  int is = 2; // hence 1..2 variables will be treated as invertible
  def A  = ratstd(I,is);
  pGBid; // polynomial form
  setring A;
  rGBid; // rational form
}
*/

/*
exParamAppelF4()
{
  // Appel F4
  LIB "ratgb.lib";
  ring r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
  matrix @D[4][4];
  @D[1,3]=1; @D[2,4]=1;
  def S=nc_algebra(1,@D);
  setring S;
  ideal I =
    x*Dx*(x*Dx+c-1) - x*(x*Dx+y*Dy+a)*(x*Dx+y*Dy+b),
    y*Dy*(y*Dy+d-1) - y*(x*Dx+y*Dy+a)*(x*Dx+y*Dy+b);
  def A = ratstd(I,2);
  pGBid; // polynomial form
  setring A;
  rGBid; // rational form
  // hence, the K(x,y) basis is {1,Dx,Dy,Dy^2}
}

// more examples:

// F1 is hard
appel F1
{
  LIB "dmodapp.lib";
  LIB "ratgb.lib";
  def A = appelF1();
  setring A;
  IAppel1;
  def F1 = ratstd(IAppel1,2);
  lead(pGBid);
  setring F1; rGBid;
}

// F2 is much easier
appel F2
{
  LIB "dmodapp.lib";
  LIB "ratgb.lib";
  def A = appelF2();
  setring A;
  IAppel2;
  def F1 = ratstd(IAppel2,2);
  lead(pGBid);
  setring F1; rGBid;
}

// F4 is feasible
appel F4
{
  LIB "dmodapp.lib";
  LIB "ratgb.lib";
  def A = appelF4();
  setring A;
  IAppel4;
  def F1 = ratstd(IAppel4,2);
  lead(pGBid);
  setring F1; rGBid;
}


*/

// Important: example for treating modules
// take two annihilators in 2 components

/*
  LIB "nctools.lib";
  ring r = (0,c),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
  // this ordering is an antiblock ordering, as it must be
  def S = Weyl(); setring S;
  // the ideal I below annihilates parametric Appel F4 function
  // where we set parameters to a=-2, b=-1 and d=0
  ideal I =
    x*Dx*(x*Dx+c-1) - x*(x*Dx+y*Dy-2)*(x*Dx+y*Dy-1),
    y*Dy*(y*Dy-1) - y*(x*Dx+y*Dy-2)*(x*Dx+y*Dy-1);

  // the ideal J below annihilates parametric Appel F4 function
  // where we set parameters to a=0, b=-1, c=0, d=0

  ideal J =
    x*Dx*(x*Dx-1) - x*(x*Dx+y*Dy)*(x*Dx+y*Dy-1),
    y*Dy*(y*Dy-1) - y*(x*Dx+y*Dy)*(x*Dx+y*Dy-1);

  module M = I*gen(1), J*gen(2);

// harder modification: M = M, Dx*gen(1) + Dy*gen(2);
// gives K(x,y)-dim 3

  int is = 2; // hence 1st and 2nd variables, that is x and y
  // will become invertible in the localization
  def A = ratstd(M,2); // main call
  pGBid; // polynomial form of the basis in the localized ring
  setring A;
  // we see from computations, that the K(x,y) dimension is 8
  rGBid; // "rational" or "localized" form of the basis
  print(matrix(kbase(rGBid)));// we see  the K(x,y) basis of the corr. module

*/