#   =================================================================================================================================  #      
                            
#  Description : Computes the border basis of a zero-dimensional ideal (which first computes a
#                Groebner basis and then converts it to a border basis in the spirit of the FGLM technique)

#   =================================================================================================================================  #

    #    Authors Contact:  
                         
          #     Amir Hashemi <amir.hashemi(AT)cc.iut.ac.ir>                 #
          #     Martin Kreuzer <Martin.Kreuzer(AT)Uni-Passau.de>            #
          #     Samira Pourkhajouei <s.pourkhajooei(AT)math.iut.ac.ir>      #

#   =================================================================================================================================  #

with(Groebner):
with(PolynomialIdeals):
#<-----------------------------------------------------------------------------
lm := proc (f, T) options operator, arrow; LeadingMonomial(f, T) end proc:
#<-----------------------------------------------------------------------------
lc := proc (f, T) options operator, arrow; LeadingCoefficient(f, T) end proc:
#<-----------------------------------------------------------------------------
ord1:= proc (a, b, T) 
local i, j, t1, t2, Aa1, Bb2;
global G;
i := a[2]; 
j := b[2]; 
t1 := a[1]; 
t2 := b[1]; 
Aa1:=expand(t1*G[i][2]);
Bb2:=expand(t2*G[j][2]);
if Aa1<>Bb2 then 
   RETURN(TestOrder(Aa1, Bb2, T)): 
else
 RETURN(evalb(i<j)): 
end if: 
end proc:
#<-----------------------------------------------------------------------------
ord2:= proc (a, b, T) 
global G:
local i, j, t1, t2, Aa1, Bb2;
i := a[2]; 
j := b[2]; 
t1 := a[1]; 
t2 := b[1]; 
Aa1:=expand(t1*G[i][2]);
Bb2:=expand(t2*G[j][2]);

 if Aa1 <> Bb2 then

            RETURN(TestOrder(Aa1, Bb2, T));

    else 
          if i=j then

            return(TestOrder(t1*G[i][5][2], t2*G[j][5][2], T));

          else

            RETURN(evalb(i <j));

          end if:

    end if;
end proc:
#<-----------------------------------------------------------------------------
Cover := proc (p1, p2, T) 
local u1, u2, v1, v2,t; 
u1 := p1[1]; 
u2 := p2[1];  
if not evalb(u1[2] = u2[2]) then 
   RETURN(false): 
elif divide(u1[1], u2[1],'q') then 
   t := q*p2[2]; 
   if TestOrder(t, p1[2],T) and t <> p1[2] then 
       RETURN(true):
   end if:
end if; 
RETURN(false):
end proc:
#<-----------------------------------------------------------------------------
regtopred := proc (p1, T)
local Flag, n, p, i, flag, u1, v1, u2, v2, t, c; 
global G:
Flag := true;
n := op(2, op(2, G)); 
p := [p1[1], p1[2]]; 
while Flag do 
   i := 1; 
   flag := false; 
   u1, v1:= p[1], p[2]; 
   while flag = false and i <= n do 
      u2, v2 := G[i][5], G[i][1]; 
      if divide(lm(v1, T), G[i][2], 'q') then 
         t := q; 
         c := lc(v1, T)/G[i][3]; 
         if ord1([t*u2[1], u2[2]], u1, T) then
            if [t*u2[1], u2[2]]<>u1 or c <> 1 then 
               p := [u1, v1-expand(c*t*v2)]; 
               flag := true; 
               if [t*u2[1], u2[2]] = u1 then 
                   p[2] := p[2]/(1-c): 
               end if:
            end if:
         end if: 
      end if; 
      i := i+1: 
   end do; 
   if flag = false then 
      Flag := false:
   end if:
end do; 
RETURN(p): 
end proc:
#<-----------------------------------------------------------------------------
firstborder:= proc (O, var)
 RETURN(`minus`({seq(seq(i*j, `in`(j, var)), `in`(i, O))}, {op(O)}));
end proc:
#<-----------------------------------------------------------------------------
gvw:= proc (F, T, var) 
local tim1, b1, NumOfTop, NumOfJP, NumOfZero, NumOfCov, n, JP, H, i, j, Lcm, t, t1T, t2T, jp, m, f, mf, flag, l, g, gs, plc, Ojp, p, tp, sg, r, flag1, jp2, t2, f2, tsig, T1, T2, JPnew, m1, m2, t1u1, t2u2, jpn, t0, i0, f0, jpairr, sp, jpair, t3, f3, tim2, b2, L, ordideal, border, BB, mm, NFB;
global G,Fideal:
tim1,b1:=kernelopts(cputime,bytesused):
NumOfTop:=0:NumOfJP:=0:
NumOfZero:=0:
NumOfCov:=0: 
n := nops(F); 
JP := [];
G:=[seq([F[i],lm(F[i],T),lc(F[i],T),i,[1,i]],i=1..nops(F))]:
G:=Array(1..n,G):   
H := [seq(seq([G[j][2], i], j = i+1 .. n), i = 1 .. n)];
for i to n do 
   for j from i+1 to n do 
      Lcm := lcm(G[i][2], G[j][2]); 
      t[i][j] := Lcm/G[i][2]; 
      t[j][i] := Lcm/G[j][2]; 
      t1T := [t[i][j]*G[i][5][1], G[i][5][2]]; 
      t2T := [t[j][i]*G[j][5][1], G[j][5][2]]; 
      if t1T<>t2T or G[i][3] <> G[j][3] then 
         if ord1(t1T, t2T, T) then 
            JP := [op(JP), [t[j][i], j]]:NumOfJP:=NumOfJP+1;
         else 
            JP := [op(JP), [t[i][j], i]]:NumOfJP:=NumOfJP+1;
         end if:
      end if:
   end do:
end do;

for jp in JP do 
    m, f := jp[1], G[jp[2]]; 
    mf := [[m*f[5][1],f[5][2]],m*f[2]]; 
    flag := false;
    l:=1: 
    while flag =false and l <= n do
         g:=G[l]: 
         gs:=[g[5],g[2]]:
         if Cover(mf, gs, T) then 
             plc := 'plc'; 
             member(jp, JP, 'plc'); 
             JP := subsop(plc = NULL, JP); 
             NumOfCov:=NumOfCov+1;
             flag := true; 
         else
             l:=l+1:
         fi:
    end do;
   
    l:=1:
    while flag = false and l<=nops(JP) do
         Ojp:=JP[l]:                               
         t,p:=Ojp[1],G[Ojp[2]]:                                                            
         tp:=[[t*p[5][1],p[5][2]],t*p[2]]:                                   
         if Cover(mf,tp,T) then                                             
             plc:='plc':                                              
             member(jp,JP,'plc'):                                                                               
             JP:=subsop(plc=NULL,JP): 
             NumOfCov:=NumOfCov+1;
             flag:=true:
         else
             l:=l+1:                                          
         fi:                             
    end do: 
end do;               

JP := sort(JP,(a,b)->ord2(a,b,T)): 

while op(JP) <> NULL do 
    t,i:=JP[1][1],JP[1][2]:                 
    p:=G[i]:                
    JP:=JP[2..-1]:                
    jp:=[[t*p[5][1],p[5][2]],expand(t*p[1]),t*p[2]]:             
    flag:=false: 
    l:=1:
    while flag=false and l<= n do 
        g:= G[l]: 
        sg:=[g[5],g[2]]:                            
        if Cover([jp[1],jp[3]],sg,T) then
            NumOfCov:=NumOfCov+1:                         
            flag:=true:                           
        else
            l:=l+1:                     
        fi:              
    end do:
  
   if  flag=false then               
       r:=regtopred(jp,T):
       NumOfTop:=NumOfTop+1:                      
       if r[2]=0 then
           NumOfZero:=NumOfZero+1:                                                   
           H:=[op(H),r[1]]:
            flag1:=false:
           l:=1:                             
           while flag1=false and l <= nops(JP) do 
               jp2:=JP[l]:                              
               t2,f2:=jp2[1],G[jp2[2]]:                                                            
               tsig:=[t2*f2[5][1],f2[5][2]]:                                   
               if tsig[2]=r[1][2] and divide(tsig[1],r[1][1],'q') then                                             
                   plc:='plc':                                              
                   member(jp2,JP,'plc'):                                                                               
                   JP:=subsop(plc=NULL,JP): 
                   NumOfCov:=NumOfCov+1:
                   flag1:=true:
               else
                   l:=l+1:                                            
               fi:                             
           end do:  
                              
       else
           r:=[r[2],lm(r[2],T),lc(r[2],T),n+1,r[1]]:  
                                                                    
           for p in G do                               
               T1:=[r[2]*p[5][1],p[5][2]]:                 
               T2:=[p[2]*r[5][1],r[5][2]]:                       
               if T1<>T2 or r[3]<>p[3] then                       
                   if ord1(T1,T2,T) then                                 
                       H:=[op(H),T2]:                          
                   else                                
                       H:=[op(H),T1]:                          
                   fi:                        
               fi:                   
           od:
           
           JPnew:=[]:                                       
           for p in G do                           
               Lcm:=lcm(p[2],r[2]):                           
               m1:=Lcm/p[2]:                          
               m2:=Lcm/r[2]:
               t1u1:=[m1*p[5][1],p[5][2]]:
               t2u2:=[m2*r[5][1],r[5][2]]:                        
               if t1u1<>t2u2 or p[3]<>r[3] then                       
                   if ord1(t1u1,t2u2,T) then                                   
                        JPnew:=[op(JPnew),[m2,n+1]]:  NumOfJP:=NumOfJP+1;                             
                   else                                    
                        JPnew:=[op(JPnew),[m1,p[4]]]: NumOfJP:=NumOfJP+1;                             
                   fi:                         
               fi:                  
           od:
                         
           n:=n+1:
           G:=Array(1..n,G): 
           G[n]:=r: 
                         
           for jpn in JPnew do
               t0,i0:=jpn[1],jpn[2]:
               f0:=G[i0]:
               jpairr:=[[t0*f0[5][1],f0[5][2]],t0*f0[2]]:
               flag1 := false;
               l:=1:
               while flag1=false and l<=n-1 do 
                    p:=G[l]: 
                    sp:=[p[5],p[2]]:
                    if Cover(jpairr, sp, T) then 
                       plc := 'plc'; 
                       member(jpn, JPnew, 'plc'); 
                       JPnew := subsop(plc = NULL, JPnew); 
                       NumOfCov:=NumOfCov+1:
                       flag1 := true;
                    else
                       l:=l+1: 
                    fi:
               end do;
              
               l:=1:
               while flag1=false and l<=nops(JP) do
                    jpair:=JP[l]:
                    t3,f3:=jpair[1],G[jpair[2]]:
                    jp2:=[[t3*f3[5][1],f3[5][2]],t3*f3[2]]:
                    if Cover(jpairr,jp2,T) then
                       plc := 'plc'; 
                       member(jpn, JPnew, 'plc'); 
                       JPnew := subsop(plc = NULL,JPnew);
                       NumOfCov:=NumOfCov+1:
                       flag1:=true:
                    else
                       l:=l+1:
                    fi:
               od:
           end do:
           JP:=[op(JP),op(JPnew)]:                          
           JP:=sort(JP,(a,b)->ord2(a,b,T)): 
       fi:                                              
   fi:      
end do:

tim2,b2:=kernelopts(cputime,bytesused):

L:=[seq(g[1],g=G)]: 
ordideal:=NormalSet(L, T)[1]:
border:=firstborder(ordideal,var):
BB:=NULL:
for mm in border do
NFB:=expand(mm-NormalForm(mm,L,T));
 BB:=BB,NFB;
od:

printf("%-1s %1s %1s %1s         : %3a %3a\n",The, cpu, time, is,tim2-tim1,(sec)):
printf("%-1s %1s %1s         : %5a %3a\n",The,used,memory,b2-b1,(bytes)):
printf("%-1s %1s %1s %1s %1s      : %5a\n",The ,number, of, zero, reduction,NumOfZero):
printf("%-1s %1s %1s %1s           : %5g\n",The ,number,of,Covercriteria,NumOfCov):
printf("%-1s %1s %1s %1s        : %5a\n",The, number, of, TopReduction,NumOfTop):
printf("%-1s %1s %1s %1s             : %5a\n",The, number,of,Jpairs,NumOfJP):

 return([BB],ordideal): 
end proc:



 

print("####################--E.1 ");
#F:=gvw([y^4 + x*y^2*z + x^2-2*x*y + y^2 + z^2, -x^3*y^2+x*y*z^3 + y^4 + x*y^2*z-2*x*y, x*y^4 + y*z^4-2*x^2*y-3],tdeg(x,y,z),[x,y,z]);

print("####################--E.2 ");
#F:=gvw([2*t^2 + u^2 + 2*x^2 + 2*y^2 + 2*z^2-u, 2*t*u + x*y + 2*t*z + 2*y*z-t, t^2 + 2*t*y + 2*u*z + 2*x*z-z, 2*t*x + 2*u*y + 2*t*z-y, 2*t + u + 2*x + 2*y + 2*z-1],tdeg(x,y,z,t,u),[x,y,z,t,u]);

print("####################--E.3 ");
F:=gvw([x+3*x*y^3+y^4+y*z^2, -x^2*z+2*y^3*z+z^2+2*y*z^2+3*x*y*z^2, 3*x^3+x*y^2+y*z^2-z*x*z^3],tdeg(x,y,z),[x,y,z]);

print("####################--E.4 ");
#F:=gvw([x^2+y^4+x^3*z+y*z-2*x*z^3, -x^2*y^2-y^3*z-z^3-3*y*z^3, y^4-x^2*z+2*y^2*z-2*x*y*z^2],tdeg(x,y,z),[x,y,z]);

print("####################--E.5 ");
#F:=gvw([2*t*x*y+y^2*z-2*x-z, t^2*x+2*t*y*z-x-2*z, 4*t*x^2*y+2*t*y^3-x^3*z+4*x*y^2*z-10*t*y+4*x^2+4*x*z-10*y^2+2, 2*t^3*y+4*t^2*x*z+4*t*y*z^2-x*z^3-10*t^2-10*t*y+4*x*z+4*z^2+2],tdeg(x,y,z,t),[x,y,z,t]);

print("####################--E.6 ");
#F:=gvw([y^4+x*y^2*z+x^2-2*x*y+y^2+z^2, y^4*x+y*z^4-2*x^2*y-3, -x^3*y^2+x*y*z^3+y^4+x*y^2*z-2*x*y],tdeg(x,y,z),[x,y,z]);

print("####################--E.7 ");
#F:=gvw([x^7-y-3*z, x*y^5-5057*z^2-2, y*z^6-x-z+14],tdeg(x,y,z),[x,y,z]);

print("####################--E.8  ");
#F:=gvw([w^5-11*w^4+41*w^3-61*w^2+30*w, 7*w^4-74*w^3+257*w^2-286*w+24*x-144, -59*w^4+614*w^3-1909*w^2+1594*w+120*y-480,107*w^4-962*w^3+2497*w^2-2002*w+120*z-120],tdeg(x,y,z,w),[x,y,z,w]);

print("####################--E.9  ");
#F:=gvw([x^4+83*x^3+73*y^2-85*z^2-427*t, y^3-x,z^3+z-t, t^3-324*z^2+94*y^2+76*x],tdeg(x,y,z,t),[x,y,z,t]);

print("####################--E.10 ");
#F:=gvw([x^2*y^2+x^2*y, y^4+2*y^3+y^2, x^3-x*y^2-x*y],tdeg(x,y),[x,y]);

print("####################--E.11 ");
#F:=gvw([x^2*y*z + x*y^2*z + x*y*z^2 + x*y*z + x*y + x*z + y*z, x^2*y^2*z + x*y^2*z^2 + x^2*y*z + x*y*z + y*z + x + z, x^2*y^2*z^2 + x^2*y^2*z + x*y^2*z + x*y*z + x*z + z + 1],tdeg(x,y,z),[x,y,z]);

print("####################--E.12 ");
#F:=gvw([x^4+83*x^3+73*y^2-85*z^2-437*t, y^3-z-t, z^3+z-t, t^4-12*z^2+77*y^2+15*x],tdeg(x,y,z,t),[x,y,z,t]);