#   =================================================================================================================================  #      
                            
 #  Description : Computes the border basis of a zero-dimensional ideal and the corresponding order ideal (Original form of Kehrein and Kreuzer algorithm)    #   

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

    #    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(PolynomialIdeals): with(Groebner):
#<-----------------------------------------------------------------------------
firstborder := proc (O, var) 
                return(`minus`({seq(seq(i*j, j in var), i in O)}, {op(O)}));
end proc:
#<-----------------------------------------------------------------------------
supp := proc (g, T)
 local N, f, A;
  N := NULL;
  f := g;
  while f <> 0 do
   A := LeadingTerm(f, T);
   f := f-A[1]*A[2];
   N := N, A[2]
  end do;
   return (N)
end proc:
#<-----------------------------------------------------------------------------
LM := proc (f, T)
 if f = 0 then
   return(0);
 else
   return(LeadingMonomial(f, T));
 end if
end proc:
#<-----------------------------------------------------------------------------  
LN := proc (f, L, T)
local m,g;
  m := LM(f, T);
  g := expand(f-m);

  return(expand(m+NormalForm(g, L, T)));

end proc:
#<-----------------------------------------------------------------------------
Finalred := proc (VV, T, O, var)
local V, VR, v, H, F, B, i;
global p;
 V := VV; VR := NULL;
 V := sort(V, proc (a, b) options operator, arrow; TestOrder(LM(a, T), LM(b, T), T) end proc);
while V <> [] do
 v := V[1];
 V := V[2 .. -1];
 H := `minus`({supp(v, T)}, {LM(v, T), op(O)});
 if nops(H) = 0 and v <> 0 then
   v:=v/LeadingCoefficient(v, T);
    VR := VR, v;
 else
   v := LN(v, [VR], T);
    if v <> 0 then
    v:=v/LeadingCoefficient(v, T);
    VR := VR, v;
    end if;
 end if;
end do;
   F := firstborder(O, var);
    B := NULL;
    for i in [VR] do
       if `subset`({LM(i, T)}, F) then
        B := B, i;
       end if;
    end do;
     return ([B]);
end proc:
#<-----------------------------------------------------------------------------
lem12 := proc (L, G, T)	
local H, n, f, i, r, V, W;
global kk,p,SAMI,pairs,stud;
H := G; V := L; W := NULL; n := 0;
pairs:=pairs+nops(H);
H := sort(H, proc (a, b) options operator, arrow; TestOrder(LM(a, T), LM(b, T), T) end proc);
while H <> [] do
  f := H[1];  H := H[2 .. -1];
  i := 1; stud:=stud+1;
   while f <> 0 and i <= nops(V) do
    if LM(f, T) = LM(V[i], T) then
     f := expand(f-LeadingCoefficient(f, T)*V[i]); 
      i := 1;
    else
      i := i+1;
    end if
   end do;
     if f <> 0 then
      f:=f/LeadingCoefficient(f, T);
       V := [op(V), f];
       W := W, f;
       else
     
       kk := kk+1;
     end if
end do;
       return ([W]);
end proc:
#<-----------------------------------------------------------------------------
posoper := proc (L, var)
  return([op(`minus`({seq(seq(expand(v*l), l = L), v = var)}, {op(L)}))]);
end proc:
#<-----------------------------------------------------------------------------
deg := proc (f, d)
  if degree(f) <= d then
    return(true);
  else
    return(false);
  end if
end proc:
#<-----------------------------------------------------------------------------
stable := proc (V1, T, var, d)
local G, W, r, L, WW, n, k, V;
 n := 1; V := V1;
 while 0 < n do
  WW := lem12(V, posoper(V, var), T);
   W := select(deg, WW, d);
   n := nops(W);
     if n <> 0 then
       V := [op(V), op(W)];
     end if;
 end do;
     return (V);
end proc:
#<-----------------------------------------------------------------------------
BorderBasis := proc (VV, T, var)
local firsttime, firstbytes, d, L, V, U, O, F, B, secondtime, secondbytes;
global kk,p,pairs,stud;
firsttime, firstbytes := kernelopts(cputime, bytesused);
 kk := 0;pairs:=0;stud:=0;
           d := max([seq(degree(i), i in VV)]);
           L := [op(randpoly(var, degree = d, dense, coeffs = proc () 1 end proc))];
           V := lem12([], VV, T);
           U := stable(V, T, var, d);
           O := `minus`({op(L)}, {seq(LM(v, T), v = U)});
           F := firstborder(O, var);
               while not `subset`(F, {op(L)}) do
                 d := d+1;
                 L := [op(randpoly(var, degree = d, dense, coeffs = proc () 1 end proc))];
                 U := stable(U, T, var, d);
                 O := `minus`({op(L)}, {seq(LM(v, T), v = U)});
                 F := firstborder(O, var);
               end do;
                 B := Finalred(U, T, O, var);
                 secondtime, secondbytes := kernelopts(cputime, bytesused);
  

    printf("%-1s %1s %1s %1s:         %3a %3a\n", The, cpu, time, is, secondtime-firsttime, sec);
    printf("%-1s %1s %1s:         %3a %3a\n", The, used, memory, secondbytes-firstbytes, bytes);
    printf("%-1s %1s %1s %1s:    %3a\n", The, number, of, Bpairs, pairs);
    printf("%-1s %1s %1s %1s %1s:    %3a\n", The, number, of, treated, pair, stud);         
    printf("%-1s %1s %1s %1s %1s:    %3a\n", The, number, of, reduction, zero, kk);

      return(B,O);

end proc:



print("####################--E.1 ");
#F:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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:=BorderBasis([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]);