source: git/Singular/LIB/swalk.lib @ 591e530

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

version="$Id$";
category="Commutative Algebra";

info="
LIBRARY: swalk.lib   Sagbi Walk Conversion Algorithm
AUTHOR: Junaid Alam Khan

PROCEDURE:
  swalk(ideal[,intvec]);   Sagbi basis of subalgebra via Sagbiwalk algorithm
";

LIB "algebra.lib";
LIB "sagbi.lib";
//////////////////////////////////////////////////////////////////////////////
proc swalk(ideal Gox, list #)
"USAGE:  swalk(i); ideal;
         swalk( i,v,w); i ideal, v,w integer vectors;
RETURN:  If list #= v,w (resp empty) then it compute the sagbi basis of
         the subalgebra defined by the generators of ideal,
         calculated via the Sagbi walk algorithm from the ordering
         (a(v),lp) (resp dp) tothe ordering (a(w),lp) ( resp lp).
EXAMPLE: example swalk; shows an example
"
{
  /* we use ring with ordering (a(...),lp,C) */
  list OSCTW    = OrderStringalp_NP("al", #);//"dp"
  string ord_str =   OSCTW[2];
  intvec icurr_weight   =   OSCTW[3]; /* original weight vector */
  intvec itarget_weight =   OSCTW[4]; /* terget weight vector */
  kill OSCTW;
  option(redSB);
  def xR = basering;
  execute("ring ostR = ("+charstr(xR)+"),("+varstr(xR)+"),"+ord_str+";");
  def new_ring = basering;
  ideal Gnew = fetch(xR, Gox);
  Gnew=sagbi(Gnew,1);
  Gnew=interreduceSd(Gnew);
  vector curr_weight=Ctv(icurr_weight);
  vector target_weight=Ctv(itarget_weight);
  ideal Gold;
  list l;
  intvec v;
  int n=0;
  while(n==0)
    {
       Gold=Gnew;
       def old_ring=new_ring;
       setring old_ring;
       number ulast;
       kill new_ring;
       if(curr_weight==target_weight){n=1;}
       else {
              l=Df(Gold);
              ulast=last(curr_weight, target_weight, l);
              vector new_weight=newvec(curr_weight,target_weight,ulast);
              vector w=cleardenom(new_weight);
              v=CT(w);
              list p= ringlist(old_ring);
              p[3][1][2]= v;
              def new_ring=ring(p);
              setring new_ring;
              ideal Gold=fetch(old_ring,Gold);
              vector curr_weight=fetch(old_ring,new_weight);
              vector target_weight=fetch(old_ring,target_weight);
              kill old_ring;
              ideal Gnew=Convert(Gold);
              Gnew=interreduceSd(Gnew);
           }
    }
   setring xR;
   ideal result = fetch(old_ring, Gnew);
   attrib(result,"isSB",1);
   return (result);
}
example
{
  "EXAMPLE:";
  //** compute a Sagbi basis of I w.r.t. lp.
  ring r = 0,(x,y), lp;
  ideal I =x2,y2,xy+y,2xy2+y3;
  intvec v=1,1;
  intvec w=1,0;
  swalk(I,v,w);
}

proc inprod(vector v,vector w)
"USAGE:  inprod(v,w); v,w vectors
RETURN:  inner product of vector v and w
EXAMPLE: example inprod; shows an example
"
{
  poly a;
  int i;
  for(i=1;i<=nvars(basering);i++)
    {
      a=a+v[i]*w[i] ;
    }
  return(a);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   vector v =[1,-1,2];
   vector w = [1,0,3];
   inprod(v,w);
}


proc df(poly f)
"USAGE:  df(f); f polynomial
RETURN:  a list of integers vectors which are the differnce of exponent vector
         of leading monomial of f with the exponent vector of
         of other monmials in f.
EXAMPLE: example df; shows an example
"
{
  list l;
  int i;
  intvec v;
  for(i=2;i<=size(f);i++)
   {
     v=leadexp(f[1])-leadexp(f[i]);
     l[size(l)+1]=v;
   }
 return(l);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   poly f = xy+z2 ;
   df(f);
}


proc Df( ideal i)
"USAGE:  Df(i); i ideal
RETURN:  a list which contain df(f), for all generators f of ideal i
EXAMPLE: example Df; shows an example
"
{
 list l;
 int j;
 for(j=1;j<=size(i);j++)
  {
   l[size(l)+1]=df(i[j]);
  }
return(l);
  }
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   ideal I = xy+z2,y3+x2y2;
   Df(I);
}


proc newvec(vector v, vector w, number u)
"USAGE:  newvec(v,w,u); v,w vectors, u number
 RETURN: the vector (1-u)v+(u)w.
 EXAMPLE: example newvec; shows an example
"
{
  vector wnew ;
  wnew=(1-u)*v+(u)*w ;
  return(wnew);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   vector v=[0,0,1];
   vector w=[1,0,0];
   number u=2/3;
   newvec(v,w,u);
}

proc CT(vector v)
"USAGE:  CT(v); v  vector
RETURN:  change the type of  vector v into integer vector.
EXAMPLE: example CT; shows an example
"
{
  intvec w ;
  int j ;
  for(j=1;j<=nvars(basering);j++)
   {
     w[j]=int(leadcoef(v[j]));
   }
  return(w);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   vector v = [2,1,3];
   CT(v);
}


proc Ctv( intvec v)
"USAGE:  Ctv(v); v integer vector
RETURN:  change the type of integer vector v into vector.
EXAMPLE: example Ctv; shows an example
"
{
vector w;
int j ;
for(j=1;j<=size(v);j++)
 {
   w=w+v[j]*gen(j);
 }
return(w);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   intvec v = 4,2,3;
   Ctv(v);
}

proc last( vector c, vector t,list l)
"USAGE: last(c,t,l); c,t vectors, l list
RETURN: a  parametric value which corresponds to vector lies along
        the path between c and t using list l of integer
        vectors. This vector is the last vector on old Sagbi cone
EXAMPLE: example last; shows an example
"
{
 number ul=1;
 int i,j,k;
 number u;
 vector v;
 for(i=1;i<=size(l);i++)
  {
    for(j=1;j<=size(l[i]);j++)
      {
        v=0;
        for(k=1;k<=size(l[i][j]);k++)
          {
            v=v+l[i][j][k]*gen(k);
          }
        poly n= inprod(c,v);
        poly q= inprod(t,v);
        number a=leadcoef(n);
        number b=leadcoef(q);
        number z=a-b;
        if(b<0)
          {
            u=a/z;
            if(u<ul) {ul=u;}
          }
        kill a,b,z,n,q ;
       }
   }
return(ul);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),(a(0,0,1),lp);
   vector v= [0,0,1];
   vector w=[1,0,0];
   ideal i=z2+xy,x2y2+y3;
    list l=Df(i);
    last(v,w,l)
}

proc initialForm(poly P)
"USAGE:  initialForm(P); P polynomial
RETURN:  sum of monomials of P with maximum w-degree
         where w is first row of matrix of a given monomial ordering
EXAMPLE: example initialForm; shows an example
"
{
 poly q;
 int i=1;
 while(deg(P[i])==deg(P[1]))
   {
     q=q+P[i];
     i++;
     if(i>size(P)) {break;}
   }
  return(q);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   poly f = x2+yz+z;
   initialForm(f);
}

proc Initial(ideal I)
"USAGE:  Initial(I); I ideal
RETURN:  an ideal which is generate by the InitialForm of the generators of I.
EXAMPLE: example Initial; shows an example
"
{
 ideal J;
 int i;
 for(i=1;i<=size(I);i++)
   {
     J[i]=initialForm(I[i]);
   }
 return(J);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),dp;
   ideal I = x+1,x+y+1;
   Initial(I);
}

proc Lift(ideal In,ideal InG,ideal Gold)
"USAGE:  Lift(In, InG, Gold); In,InG, Gold ideals
         Gold given by Sagbi basis {g_1,...,g_t}
         In given by tne initial forms In(g_1),...,In(g_t)
         InG={h_1,...,h_s} a Sagbi basis of In
RETURN:  P_j, a polynomial in K[y_1,..,y_t] such that
         h_j=P_j(In(g_1),...,In_(g_t)) and return f_j=P_j(g_1,...,g_t)
EXAMPLE: example Lift; shows an example
"
{
  int i;
  ideal J;
  for(i=1;i<=size(InG);i++)
  {
    poly f=InG[i];
    list l=algebra_containment(f,In,1);
    def s=l[2];
    map F=s,maxideal(1),Gold ;
    poly g=F(check);
    ideal k=g;
    J=J+k;
    kill g,l,s,F,f,k;
  }
  return(J);
 }
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),(a(2,0,3),lp);
   ideal In = xy+z2,x2y2;
   ideal InG=xy+z2,x2y2,xyz2+1/2z4;
   ideal Gold=xy+z2,y3+x2y2;
   Lift(In,InG,Gold);
}


proc Convert( ideal Gold )
"USAGE: Convert(I); I ideal
RETURN: Convert old Sagbi basis into new Sagbi basis
EXAMPLE: example Convert; shows an example
"
{
 ideal In=Initial(Gold);
 ideal InG=sagbi(In,1)+In;
 ideal Gnew=Lift(In,InG,Gold);
 return(Gnew);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   ideal I=xy+z2, y3+x2y2;
   Convert(I);
}

proc interreduceSd(ideal I)
"USAGE:  interreduceSd(I); I ideal
RETURN:  interreduceSd the set of generators if I with
         respect to a given term ordering
EXAMPLE: example interreduceSd; shows an example
"
{
  list l,M;
  ideal J,B;
  int i,j,k;
  poly f;
  for(k=1;k<=size(I);k++)
    {l[k]=I[k];}
  for(j=1;j<=size(l);j++)
   {
     f=l[j];
     M=delete(l,j);
     for(i=1;i<=size(M);i++)
       { B[i]=M[i];}
     f=sagbiNF(f,B,1);
     J=J+f;
   }
  return(J);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,z),lp;
   ideal I = xy+z2,x2y2+y3;
   interreduceSd(I);
}

static proc OrderStringalp(string Wpal,list #)
{
  int n= nvars(basering);
  string order_str;
  intvec curr_weight, target_weight;
  curr_weight = system("Mivdp",n);
  target_weight = system("Mivlp",n);

   if(size(#) != 0)
   {
     if(size(#) == 1)
     {
       if(typeof(#[1]) == "intvec")
       {
         if(Wpal == "al"){
           order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
         }
         else {
           order_str = "(Wp("+string(#[1])+"),C)";
         }
         curr_weight = #[1];
       }
       else
       {
        if(typeof(#[1]) == "string")
        {
          if(#[1] == "Dp") {
            order_str = "Dp";
          }
          else {
            order_str = "dp";
          }
        }
        else {
          order_str = "dp";
        }
     }
    }
    else
    {
     if(size(#) == 2)
     {
       if(typeof(#[2]) == "intvec")
       {
         target_weight = #[2];
       }
       if(typeof(#[1]) == "intvec")
       {
         if(Wpal == "al"){
           order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
         }
         else {
           order_str = "(Wp("+string(#[1])+"),C)";
         }
         curr_weight = #[1];
       }
       else
       {
        if(typeof(#[1]) == "string")
        {
          if(#[1] == "Dp") {
            order_str = "Dp";
           }
           else {
              order_str = "dp";
           }
        }
        else {
           order_str = "dp";
        }
      }
     }
    }
   }
   else {
     order_str = "dp";
   }
   list result;
   result[1] = order_str;
   result[2] = curr_weight;
   result[3] = target_weight;
   return(result);
}

static proc OrderStringalp_NP(string Wpal,list #)
{
  int n= nvars(basering);
  string order_str = "dp";
  int nP = 1;// call LatsGB to compute the wanted GB  by pwalk
  intvec curr_weight = system("Mivdp",n); //define (1,1,...,1)
  intvec target_weight = system("Mivlp",n); //define (1,0,...,0)
  if(size(#) != 0)
  {
    if(size(#) == 1)
    {
      if(typeof(#[1]) == "intvec") {
        curr_weight = #[1];

        if(Wpal == "al"){
          order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
        }
        else {
          order_str = "(Wp("+string(#[1])+"),C)";
        }
      }
      else {
        if(typeof(#[1]) == "int"){
          nP = #[1];
        }
        else {
          print("// ** the input must be \"(ideal, int)\" or ");
          print("// **                   \"(ideal, intvec)\"");
          print("// ** a lex. GB will be computed from \"dp\" to \"lp\"");
        }
      }
    }
    else
    {
     if(size(#) == 2)
     {
       if(typeof(#[1]) == "intvec" and typeof(#[2]) == "int"){
         curr_weight = #[1];

         if(Wpal == "al"){
           order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
         }
         else {
           order_str = "(Wp("+string(#[1])+"),C)";
         }
       }
       else{
         if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec"){
           curr_weight = #[1];
           target_weight = #[2];

           if(Wpal == "al"){
             order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
           }
           else {
             order_str = "(Wp("+string(#[1])+"),C)";
           }
         }
         else{
           print("// ** the input  must be \"(ideal,intvec,int)\" or ");
           print("// **                    \"(ideal,intvec,intvec)\"");
           print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");
         }
       }
     }
     else {
       if(size(#) == 3) {
         if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec" and
            typeof(#[3]) == "int")
         {
           curr_weight = #[1];
           target_weight = #[2];
           nP = #[3];
           if(Wpal == "al"){
             order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)";
           }
           else {
             order_str = "(Wp("+string(#[1])+"),C)";
           }
         }
         else{
           print("// ** the input must be \"(ideal,intvec,intvec,int)\"");
           print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");

         }
       }
       else{
         print("// ** The given input is wrong");
         print("// ** and a lex. GB will be computed from \"dp\" to \"lp\"");
       }
     }
    }
  }
  list result;
  result[1] = nP;
  result[2] = order_str;
  result[3] = curr_weight;
  result[4] = target_weight;
  return(result);
}