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

THEORY: Let A be an operator algebra with R = K[x1,.,xN] as subring.
The operators are usually denoted by {d1,..,dM}.
@* Assume, that A is a G-algebra, then the set 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 sa=bt.
@* In other words, one can transform any left fraction into the right fraction. The algebra 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.

PROCEDURES:
ratstd(ideal I, int n);   compute Groebner basis and dimensions in Ore localization of the basering w.r.t. first n variables

SUPPORT: SpezialForschungsBereich F1301 of the Austrian FWF and
@* Transnational Access Program of RISC Linz, Austria
"

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 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 rGBid in C
represents the rational form of the output ideal 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,(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
}

/*
proc 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;
}


*/