//////////////////////////////////////////////////////////////////////////////
version="$Id: ratgb.lib,v 1.12 2008-02-26 23:36:12 levandov Exp $";
category="Noncommutative";
info="
LIBRARY: ratgb.lib  Groebner bases in Ore localizations
AUTHOR: Viktor Levandovskyy,     levandov@risc.uni-linz.ac.at

PROCEDURES:
ratstd(ideal I, int n);   compute Groebner basis in Ore localization of the basering wrt first n variables

SUPPORT: SpezialForschungsBereich F1301 of the Austrian FWF
"

LIB "poly.lib";

//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 ordering
@*     for the variables in the second block
NOTE: the output ring O is commutative. The ideal rGBid in O
represents the rational form of the output ideal pGBid in the basering.
@* During the computation, the D-dimension of I and the corresponding
vector space D-dimension of I are computed and printed out.
@* Setting optional integer eng to 1, slimgb is taken as Groebner engine
TRACING: 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;

  // assume the ordering is the block with (a(0:is),ORD)
  // set up ORD as the ordering
  //  L; "RL:"; RL;
  if (size(L[3]) != 3)
  {
    "note: strange ordering";
  }
  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);

  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) wrt. 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,1),dp);
  matrix D[4][4];
  D[1,3] = K;
  D[2,4] = N;
  def S = nc_algebra(1,D);
  setring S;
  ideal I = (k+1)*K - (n-k), (n-k+1)*N - (n+1);
  int is = 1;
  def A  = ratstd(I,is);
  pGBid; // polynomial form
  setring A;
  rGBid; // rational form
}

static proc exParam()
{
  // Appel F4
  LIB "ratgb.lib";
  ring r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),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
}