LIB "matrix.lib";

ring r=0,(x,y),lp;

ideal i=x4+x2+xy3+2,x2+y2-1;

ideal g=std(i);
//option(redSB);
g;
dim(g);
vdim(g);

ideal B=kbase(g);
int   L=size(B);


matrix m[L][1];
matrix column[L][1];

int n=0;          //power of the ring variable for which we want to calculate the elimination ideal
int flag=1;       //flag is set to zero, when the reduced powers become linearly dependent
while(flag==1)
{

poly f=reduce(x^n,g);

for(int counter=1;counter<=L;counter=counter+1)
{
  column[counter,1]=jet(f/B[counter],0);
}

if (n==0)
    {m=column;}
else
     {m=concat(m,column);}

if (rank(m)<(n+1))
      {flag=0;}

n=n+1;

}

print(m);

///////////////////////////////////////////////////
option(noredefine);

int nc=ncols(m);

ring r2=0,(x(1..nc)),lp;    //number of variables must be the same as the number of columns of the matrix for which we want the reduced row echelon form

map f=r,x(1),x(2);
matrix m2=f(m);


print("m2=");print(m2);

int nr=nrows(m2);
int nc=ncols(m2);
int j;
int k;
for(j=1;j<=nr;j=j+1){for(k=1;k<=nc;k=k+1){m2[j,k]=m2[j,k]*x(k);}};
kill j;
kill k;
//print(m2);

matrix ONES[nc][1];
for(int j=1;j<=nc;j=j+1){ONES[j,1]=1;};
kill j;

ideal i=m2*ONES;
option(redSB);
ideal g=std(i);
option(redSB);
int L=size(g);
for(int j=1;j<=L;j=j+1) {g[j]=g[j]/leadcoef(g[j]);}
kill j;

matrix Output[L][nc];

int k;
for (int j=L;j>0;j=j-1){for (k=1;k<=nc;k=k+1){Output[j,k]=jet(g[L+1-j]/x(k),0);}};
kill j;
kill k;


print("Output=");print(Output);print("Output is the reduced row echelon form of matrix m");

setring r;

map phi=r2,x;
matrix rref=phi(Output);

print(rref);

poly EliminationIdeal=x^(nc-1);

for (int n=1;n<=(nc-1);n=n+1)
    {EliminationIdeal=EliminationIdeal+(-rref[n,nc])*x^(n-1);};
print(EliminationIdeal);


//ring r=0,(x,y),dp; ideal i=x3-yx+1,x2+y2-1;---->x6+x4+2x3-x2+1
//ring r=0,(x,y),dp; ideal i=x3-y2+xy,x+y-1; ---->x3-2x2+3x-1
//ring r=0,(x,y),lp; ideal i=x4+x2+xy3+2,x2+y2-1;---->x8-1/2x6+4x4+3/2x2+2