rstandard.lib

///////////////////////////////////////////////////////////////////////////////
version="version standard.lib 4.1.2.0 Feb_2019 ";
category="Commutative Algebra";
info="
LIBRARY: rstandard.lib     Computes Janet bases and border bases for ideals

AUTHORS: Shamsa Kanwal           lotus_zone16@yahoo.com
         Gerhard Pfister         pfister@mathematik.uni-kl.de

OVERVIEW: Computing Janet bases and border bases for any ordering
          using the idea of r-standard bases (defined by V. Gerdt)

REFERENCES:
   [1] A. Kehrein, M. Kreuzer, L. Robbiano: An algebrists view on border bases, in:
       A. Dickenstein and I. Emiris (eds.), Solving Polynomial Equations:
       Foundations, Algorithms and Applications, Springer, Heidelberg 2005, 169-202.

   [2] V.P. Gerdt: Involute Algorithms for Computing Groebner Bases, In
       Computational Commutative and Non-Computational Algebra Geometry,
       S.Conjocaru, G. Pfister and V. Ufnarovski (Eds.), NATO Science Series,105
       Press 2005, 199-255.


PROCEDURES:
   borderBasis(I);  computes a border basis of the ideal I
   modBorder(I);    computes a border basis of the ideal I using modular methods
   rJanet(I);       computes a Janet basis of the ideal I
   modJanet(I);     computes a Janet basis of the ideal I using modular methods

KEYWORDS: border basis; janet basis
";
LIB "general.lib";
LIB "modstd.lib";

//=================================  procedures  ==============================
//================================= border bases ==============================


static proc orderIdeal(ideal I)
{
//returns the order ideal, the corners and the border of the
//order ideal
   I=std(I);
   int d=vdim(I);
   if(d==-1){ERROR(" ideal is not zero dimensional ");}
   ideal O=kbase(I);   //this is the order ideal
   ideal C=lead(I);    //this are the corners
   ideal B=border(C);  //this is the border of the order ideal
   list L=O,C,B,I;
   return(L);
}

static proc border(ideal L)
{
//L a monomial ideal
//returns the border of the order ideal corresponding to L
   ideal O=kbase(L);
   ideal Omax=maxO(O);
   ideal B=L;
   poly p;
   list C;
   int i,j,k;
   for(i=1;i<=size(L);i++)
   {
      for(j=1;j<=size(Omax);j++)
      {
         for(k=1;k<=nvars(basering);k++)
         {
            p=(var(k)*Omax[j])/L[i];
            if(p!=0){break;}
         }
         if(k<=nvars(basering))
         {
            C=coneM(p);
            for(k=1;k<=size(C);k++)
            {
               B[size(B)+1]=C[k]*L[i];
            }
         }
      }
   }
   B=simplify(B,4);
   B=simplify(B,2);
   B=sort(B)[1];
   return(B);
}

static proc coneM(poly p)
{
//computes all monomials different from 1 dividing p
  list L=p;
  int i,j,k;
  poly q;
  k=1;
  ideal I;
  while(k)
  {
     k=0;
     i++;
     I=maxideal(i);
     for(j=1;j<=size(I);j++)
     {
       q=p/I[j];
       if((q!=0)&&(q!=1))
       {
          k=1;
          L[size(L)+1]=q;
       }
     }
  }
  return(L);
}

static proc maxO(ideal O)
{
//O is a monomial ideal
//computes the maximal elements in O with respect to division
    def R=basering;
    def S=changeord(list(list("ds",nvars(R))));
    setring S;
    ideal O=imap(R,O);
    O=sort(O)[1];
    setring R;
    O=imap(S,O);
    int i=1;
    int j;
    int n=size(O);
    while(i<ncols(O))
    {
       j=i+1;
       while(j<=ncols(O))
       {
          if(O[j]!=0)
          {
             if(O[i]/O[j]!=0)
             {
                O[j]=0;
             }
          }
          j++;
       }
       O=simplify(O,2);
       i++;
    }
    return(O);
}

/////////////////////////////////////////////////////////////////////////////
proc borderBasis(ideal I)
"USAGE:   borderBasis(I); I is an ideal.
RETURN:   ideal, a border basis for I.
PURPOSE:  Computes a border basis for the ideal given by the generators in I.
SEE ALSO: modBorder
KEYWORDS: order ideal, border basis
EXAMPLE:  example borderBasis; shows an example
"
{
   list L=orderIdeal(I);
   I=L[4];
   ideal B=L[3];
   poly P;
   int i;
   if(ringlist(basering)[3][1][1]=="ds"||ringlist(basering)[3][1][1]=="Ds")
   {
      for(i=1;i<=size(B);i++)
      {
         P=system("reduce_bound",B[i],I,deg(highcorner(I))+1);
         B[i]=B[i]-P;
      }
   }
   else
   {
      if(ord_test(basering))
      {
         for(i=1;i<=size(B);i++)
         {
            P=reduce(B[i],I);
            B[i]=B[i]-P;
         }
      }
      else
      {
         L=division(B[i],I);
         B[i]=L[3][1,1]*B[i]-L[2][1];
      }
   }
   attrib(B,"isSB",1);
   return(B);
}
example
{ "EXAMPLE:";
   ring R=32003,(x,y,z),ds;
   poly f=x3y+x5+x3y2+2x2y3+x2yz2+xy5+x12+y16+z20;
   ideal i= jacob(f);
   i=i,f;
   ideal j=borderBasis(i);  j;
}

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

proc modBorder(ideal I, list #)
"USAGE:   modBorder(I,i); I is an ideal, i an integer.
RETURN:   ideal, a border basis for I using modular methods.
PURPOSE:  Computes a border basis for the ideal given by the generators in I
          using modular techniques.
          If second argument is 0 then the result is not verified.
SEE ALSO: borderBasis
KEYWORDS: order ideal, border basis
EXAMPLE:  example modBorder; shows an example
"
{
    int exactness = 1;
    if (size(#) > 0) {
        exactness = #[1];
    }
    /* save options */
    intvec opt = option(get);
    option(redSB);
    /* choose the right command */
    string command = "borderBasis";
    /* call modular() */
    if (exactness) {
        I = modular(command, list(I), Modstd::primeTest_std,
            Modstd::deleteUnluckyPrimes_std, Modstd::pTest_std,finalTest_std_new);
    }
    else {
        I = modular(command, list(I), Modstd::primeTest_std,
            Modstd::deleteUnluckyPrimes_std,Modstd::pTest_std);
    }
    /* return the result */
    attrib(I, "isSB", 1);
    option(set, opt);
    return(I);
}
example
{ "EXAMPLE:";
   ring R=0,(x,y,z),ds;
   poly f=x3y+x5+x3y2+2x2y3+x2yz2+xy5+x12+y16+z20;
   ideal i= jacob(f);
   i=i,f;
   ideal j=modBorder(i,1);  j;

   ring S=0,(x,y,z),ds;
   ideal i=
   3x2+6xy+3y2+6xz+6yz+3z2+3x2yz+4xy2z+y3z+4xyz2+2y2z2+yz3+10x9,
   3x2+6xy+3y2+6xz+6yz+3z2+x3z+4x2yz+3xy2z+2x2z2+4xyz2+xz3+10y9,
   3x2+6xy+3y2+6xz+6yz+3z2+x3y+2x2y2+xy3+4x2yz+4xy2z+3xyz2+10z9;
   ideal j=modBorder(i,0);  j;
}

//=================================== Janet bases ==============================


/////////////////////////////////////////////////////////////////////////////
proc rJanet(ideal I)
"USAGE:   rJanet(I); I is an ideal.
RETURN:   ideal, a Janet basis for I.
PURPOSE:  Computes a Janet basis for the ideal given by the generators in I for any ordering.
SEE ALSO: modJanet
KEYWORDS: Janet division, Janet basis
EXAMPLE:  example rJanet; shows an example
"
{
   def R=basering;
   int i,d;
   poly p;

   I=std(I);
   d=dim(I);
   if(I[1]==1){return(ideal(1));}
   ideal J=lead(I);
   intvec v;
   v[nvars(R)]=0;
   v=v+1;
   list rl=ringlist(R);
   rl[3]=list(list("dp",v));
   ring S=ring(rl);
   ideal J=imap(R,J);
   J=janet(J);
   J=simplify(J,1+2+4+8+16);
   setring R;
   J=imap(S,J);
   list L;
   poly P;
   if((ringlist(basering)[3][1][1]=="ds"||ringlist(basering)[3][1][1]=="Ds")&&(d==0))
   {
      for(i=1;i<=size(J);i++)
      {
         P=system("reduce_bound",J[i],I,deg(highcorner(I))+1);
         J[i]=J[i]-P;
      }
   }
   else
   {
      if(ord_test(basering))
      {
         for(i=1;i<=size(J);i++)
         {
            P=reduce(J[i],I);
            J[i]=J[i]-P;
         }
      }
      else
      {
         L=division(J[i],I);
         J[i]=L[3][1,1]*J[i]-L[2][1];
      }
   }
   attrib(J,"isSB",1);
   return(J);
}
example
{ "EXAMPLE:";
   ring R= 32003,(a,b,c,d,e,f,g),dp;
   ideal i=
   a+b+c+d+e+f+g,
   ab+bc+cd+de+ef+fg+ga,
   abc+bcd+cde+fde+efg+fga+gab,
   abcd+bcde+cdef+defg+efga+fgab+gabc,
   abcde+bcdef+cdefg+defga+efgab+fgabc+gabcd,
   abcdef+bcdefg+cdefga+defgab+efgabc+fgabcd+gabcde,
   abcdefg-1;
   ideal j = rJanet(i);   j;
}

/////////////////////////////////////////////////////////////////////////////
proc modJanet(ideal I, list #)
"USAGE:   modJanet(I,i); I is an ideal, i an integer (optional).
RETURN:   ideal, a Janet basis for I using modular methods.
PURPOSE:  Computes a Janet basis for the ideal given by the generators in I
          using modular techniques.
          If second argument is 0 then the result is not verified.
SEE ALSO: rJanet
KEYWORDS: Janet basis
EXAMPLE:  example Janet; shows an example
"
{
    int exactness = 1;
    if (size(#) > 0) {
        exactness = #[1];
    }
    /* save options */
    intvec opt = option(get);
    option(redSB);
    /* choose the right command */
    string command = "rJanet";
    /* call modular() */
    if (exactness) {
        I = modular(command, list(I), Modstd::primeTest_std,
            Modstd::deleteUnluckyPrimes_std, Modstd::pTest_std,finalTest_std_new);
    }
    else {
        I = modular(command, list(I), Modstd::primeTest_std,
            Modstd::deleteUnluckyPrimes_std,Modstd::pTest_std);
    }
    /* return the result */
    attrib(I, "isSB", 1);
    option(set, opt);
    return(I);
}
example
{ "EXAMPLE:";
   ring R=0,(t,x,y,z),ds;
   ideal i=
   5t3x2z+2t2y3x5,
   7y+4x2y+y2x+2zt,
   3tz+3yz2+2yz4;
   ideal j=modJanet(i);  j;

   ring S=0,(x,y,z),dp;
   poly p1 =x2y*(47x5y7z3+28xy5z8+63+91x5y3z7);
   poly p2 =xyz*(57y6+21x2yz9+51y2z2+15x2z4);
   poly p3 =xy4z*(74y+32x6z7+53x5y2z+17x2y3z);
   poly p4 =y3z*(21x2z6+32x10y6z5+23x5y5z7+27y2);
   poly p5 =xz*(36y2z2+81x9y10+19x2y5z4+79x4z6);
   ideal i =p1,p2,p3,p4,p5;
   ideal j=modJanet(i,0);  j;
}
/////////////////////////////////////////////////////////////////////////////
static proc minimizeS(ideal I)
{
   int i,j;
   ideal K;
   ideal J=std(lead(I));
   I=sort(I)[1];
   J=sort(J)[1];
   while(i<size(J))
   {
      i++;
      while(j<size(I))
      {
         j++;
         if(leadmonom(I[j])==J[i])
         {
            K[size(K)+1]=I[j];
            break;
         }
      }
   }
   return(K);
}

static proc finalTest_std_new(string command, alias list args, ideal result)
{
   result=minimizeS(result);
   attrib(result,"isSB",1);
   return(Modstd::finalTest_std(command, args,result));
}

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

/*
//=====  Examples (http://symbolicdata.org/XMLResources/IntPS/): ============

//Cyclic_7.xml
ring r18= 32003,(a,b,c,d,e,f,g),dp;
ideal i=
a+b+c+d+e+f+g,
ab+bc+cd+de+ef+fg+ga,
abc+bcd+cde+fde+efg+fga+gab,
abcd+bcde+cdef+defg+efga+fgab+gabc,
abcde+bcdef+cdefg+defga+efgab+fgabc+gabcd,
abcdef+bcdefg+cdefga+defgab+efgabc+fgabcd+gabcde,
abcdefg-1;

// Milnor1.xml
ring s8=32003,(x,y,z),dp;
int a =60;
int b =40;
int c =20;
int t =1;
poly f= x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2;
ideal i= jacob(f);

// Tjurina1.xml
ring s10=32003,(x,y,z),dp;
int a=30;
poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
ideal i= jacob(f);
i=i,f;

// random1.xml
ring s15=32003,(x,y,z),dp;
poly p1 =x2y*(47x5y7z3+28xy5z8+63+91x5y3z7);
poly p2 =xyz*(57y6+21x2yz9+51y2z2+15x2z4);
poly p3 =xy4z*(74y+32x6z7+53x5y2z+17x2y3z);
poly p4 =y3z*(21x2z6+32x10y6z5+23x5y5z7+27y2);
poly p5 =xz*(36y2z2+81x9y10+19x2y5z4+79x4z6);
ideal i =p1,p2,p3,p4,p5;

// random2.xml
ring s16=32003,(x,y,z),dp;
poly p1 =xy2z*(47x6y6z2+28y5+63y3z3+91xyz6+57z2);
poly p2 =x3y2z*(21y4z7+51z9+15x4y2z7+74x7y7z5+32x7y8);
poly p3 =xy*(53x+17y2z3+21x5y5z2+32x4y7z8+23x2y4z5);
poly p4 =z*(27x7y4z6+36x3yz4+81x3y7z+19y5+79x5z6);
ideal i =p1,p2,p3,p4;

// Singular.schwarz_6.xml (dehomogenized: h=1)
ring r16= 32003,(x(1..6)),lp;
poly s1= -1*x(1)^2-1*x(1)-2*x(2)*x(6)-2*x(3)*x(5)-1*x(4)^2;
poly s2= -2*x(1)*x(2)-1*x(2)-2*x(3)*x(6)-2*x(4)*x(5);
poly s3= -2*x(1)*x(3)-1*x(2)^2-1*x(3)-2*x(4)*x(6)-1*x(5)^2;
poly s4= -2*x(1)*x(4)-2*x(2)*x(3)-1*x(4)-2*x(5)*x(6);
poly s5= -2*x(1)*x(5)-2*x(2)*x(4)-1*x(3)^2-1*x(5)-1*x(6)^2;
poly s6= -2*x(1)*x(6)-2*x(2)*x(5)-2*x(3)*x(4)-1*x(6);
ideal i =s1,s2,s3,s4,s5,s6;

//Twomat3.xml
ring r1=32003,(a,b,c,d,e,f,g,h,i,A,B,C,D,E,F,G,H,I),dp;
ideal i=
-dB-gC+bD+cG,
-bA+aB-eB-hC+bE+cH,
-cA-fB+aC-iC+bF+cI,
dA-aD+eD-dE-gF+fG,
dB-bD-hF+fH,
dC-cD-fE+eF-iF+fI,
gA+hD-aG+iG-dH-gI,
gB+hE-bG-eH+iH-hIgC+hF-cG-fH;

// Klein.xml
ring r12 = 32003,(x,y,z,t,u),dp;
ideal i= x6+y6+522*(x^5*y-x*y^5)-10005*(x^4*y^2+x^2*y^4)-z,
-(x^4+y^4)+228*(x^3*y-x*y^3)-494*x^2*y^2-t,
x*y*(x*2+11*x*y-y^2)*5-u;

// Meintjes.xml
ring r15 = 32003,(x,y,z,t,u,a,b,c,d,e,f,g),dp;
ideal i=
xy+x-3u,
yz2+2y2a+ytb+yzd+2xy+yc-ug+x,
2yz2+yzd+2z2f+ze-8u,
ytb+2t2+4ug,
yz2+y2a+ytb+yzd+z2f+xy+t2+yc+ze+x-1;

// Wilfred.xml
ring w7=32003,(x(1..16)),dp;
ideal i=
x(1)+x(2)+x(3)+x(4),
x(1)^2*x(2)^2+x(2)*x(3)*x(5)*x(6),
x(3)*x(4)*x(9)*x(10)+x(4)*x(5)*x(13)*x(14),
x(1)*x(2)^2*x(3)+x(2)*x(3)*x(6)*x(7),
x(3)*x(4)*x(11)*x(12)+x(4)*x(5)*x(15)*x(16),
x(7)*x(8)*x(9)*x(10)+x(8)*x(9)*x(13)*x(14),
x(7)*x(8)*x(11)*x(12)+x(8)*x(9)*x(15)*x(16),
x(1)*x(2)*x(9)*x(10)+x(5)*x(6)*x(10)*x(11),
x(2)*x(3)*x(9)*x(10)+x(6)*x(7)*x(10)*x(11),
x(11)*x(12)*x(15)*x(16)+x(1)*x(15)*x(16)^2,
x(4)*x(5)*x(13)*x(14)+x(8)*x(9)*x(14)*x(15)+x(12)*x(13)*x(15)*x(16)+
x(1)^2*x(16)^2;

// Behnke.xml
ring bn=32003,(a,b,c,d,e),dp;
int n=7;
ideal i=
a^n-b^n,
b^n-c^n,
c^n-d^n,
d^n-e^n,
a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;

//Paris.ilias13.xml
ring r=32003,(S1,s1,d1,S2,D2,s2,d2),dp;
ideal i=-8*d1*D2+2*S2*s2-2*D2*d2-8*S2-4*s2+16,
8*d1*D2-2*S2*s2+2*D2*d2+8*S2+4*s2-16,
-8*d1*S2-16*s1*D2+2*D2*s2-2*S2*d2+16*d1+8*D2+4*d2,
8*d1*S2+16*s1*D2-2*D2*s2+2*S2*d2-16*d1-8*D2-4*d2,
-8*S1*S2+4*S2^2-20*D2^2+16*S1-16,
8*S1*S2-4*S2^2+20*D2^2-16*S1+16,
-16*d1^2-8*s1*s2+s2^2-8*d1*d2-d2^2+32*s1-16,
16*d1^2+8*s1*s2-s2^2+8*d1*d2+d2^2-32*s1+16,
-16*S1*d1+8*d1*S2-10*D2*s2-4*S1*d2+2*S2*d2+16*d1+40*D2+4*d2,
16*S1*d1-8*d1*S2+10*D2*s2+4*S1*d2-2*S2*d2-16*d1-40*D2-4*d2,
-32*S1*s1+16*s1*S2+40*d1*D2+4*S1*s2-2*S2*s2+10*D2*d2+16*S1+32*s1-8*S2-4*s2-16,
32*S1*s1-16*s1*S2-40*d1*D2-4*S1*s2+2*S2*s2-10*D2*d2-16*S1-32*s1+8*S2+4*s2+16,
S2^2*s2^3-2*S2*D2*s2^3+D2^2*s2^3+3*S2^2*s2^2*d2-6*S2*D2*s2^2*d2+3*D2^2*s2^2*d2
  +3*S2^2*s2*d2^2-6*S2*D2*s2*d2^2+3*D2^2*s2*d2^2+S2^2*d2^3-2*S2*D2*d2^3
  +D2^2*d2^3-32,
S2^2*s2^3+2*S2*D2*s2^3+D2^2*s2^3-3*S2^2*s2^2*d2-6*S2*D2*s2^2*d2-3*D2^2*s2^2*d2
  +3*S2^2*s2*d2^2+6*S2*D2*s2*d2^2+3*D2^2*s2*d2^2-S2^2*d2^3-2*S2*D2*d2^3
  -D2^2*d2^3-32,
S1^2*s1^3-3*S1^2*s1^2*d1+3*S1^2*s1*d1^2-S1^2*d1^3+4*S1*s1^3*D2-12*S1*s1^2*d1*D2
  +12*S1*s1*d1^2*D2-4*S1*d1^3*D2+4*s1^3*D2^2-12*s1^2*d1*D2^2+12*s1*d1^2*D2^2
  -4*d1^3*D2^2-32,
S1^2*s1^3+3*S1^2*s1^2*d1+3*S1^2*s1*d1^2+S1^2*d1^3-4*S1*s1^3*D2-12*S1*s1^2*d1*D2
  -12*S1*s1*d1^2*D2-4*S1*d1^3*D2+4*s1^3*D2^2+12*s1^2*d1*D2^2+12*s1*d1^2*D2^2
  +4*d1^3*D2^2-32;


// Pfister_2.xlm
ring r=32003,(x2,x3,x4,x5,y1,y2,y3,y4,z1,z2,z3,z4,z5),dp;
ideal i=x2^3+x3^2+x4^2+x5^2+y1^2+y2^2+y3^2+y4^2+z1^2+z2^2+z3^2+z4^2+z5^2,
y1^2+z1^2-1,x2^2+y2^2+z2^2-1,x3^2+y3^2+z3^2-1,x4^2+y4^2+z4^2-1,
x5^2+z5^2-1,y1*y2+z1*z2,x2*x3+y2*y3+z2*z3,x3*x4+y3*y4+z3*z4,x4*x5+z4*z5,
x2+x3+x4+x5+1,y1+y2+y3+y4-1,z1+z2+z3+z4+z5;

// Singular.gerhard_3.xml
ring r=32003,(w,x,y,z),dp;
ideal i=3*w^5*x^5*y^5*z^2+18*z^17,
5*w^10*x^8*y^4+6*w^11*x^6*y^5+16*w^10*x^5*y^7+10*w^9*x^4*y^9+23*y^22+5*w^10*x^5*y^4*z^3,
24*x^23+8*w^11*x^7*y^5+6*w^12*x^5*y^6+10*w^11*x^4*y^8+4*w^10*x^3*y^10+5*w^11*x^4*y^5*z^3;

//Singular.gerhard_1.xml  (dehomogenized w=1)
ring r1=32003,(t,x,y,z),ds;
ideal i=
5t3x2z+2t2y3x5,
7y+4x2y+y2x+2zt,
3tz+3yz2+2yz4;

// Milnor3.xml
ring r3=32003,(x,y,z),ds;
int a =33;
int b =30;
int c =10;
int t =1;
poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2;
ideal i= jacob(f);

// Tjurina1.xml
ring r4=0,(x,y,z),ds;
int a=30;
poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
ideal i= jacob(f);

// Milnor4.xml
ring r6=32003,(x,y,z),ds;
int a =12;
int b =16;
int c =20;
int alpha =5;
int beta= 5;
int t= 1;
poly f =x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)
+x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
ideal i= jacob(f);
i=i,f;

// Steidel_1.xml
ring r=0,(x,y,z),ds;
ideal i=3*x^2+10*x^9+6*x*y+3*y^2+6*x*z+6*y*z+3*x^2*y*z+4*x*y^2*z+y^3*z
  +3*z^2+4*x*y*z^2+2*y^2*z^2+y*z^3,
3*x^2+6*x*y+3*y^2+10*y^9+6*x*z+x^3*z+6*y*z+4*x^2*y*z+3*x*y^2*z+3*z^2
  +2*x^2*z^2+4*x*y*z^2+x*z^3,
3*x^2+6*x*y+x^3*y+3*y^2+2*x^2*y^2+x*y^3+6*x*z+6*y*z+4*x^2*y*z
  +4*x*y^2*z+3*z^2+3*x*y*z^2+10*z^9;

*/