with(Algebraic):
with(Groebner):
with(PolynomialIdeals):
with(linalg):
#*************************************This is the code of Wang-Sun algorithm for factoring univariate polynomial over algebraic fields
#*************************************This procedure returns true if the input is squrefree and false, otherwise.
issqfr:=proc(f)
local V,m,n,L,i:
    V:=IdealInfo[Variables](<f>):
    L:=[seq(p[2] , p in sqrfree(f,[y])[2])]:
    n:=nops(L):m:=0:
    for i from 1 to n while L[i]<2 do
        m:=m+1:     
    od:
    RETURN(evalb(m=n)):
end:
#***********************************This is the standard algorithm.
stan:=proc(f,J,vars)
local j,fff,ggg,hhh,P,p,G0,V,tord,ns,rv,flag,sq,M,chp,lambda,N,G,pol,aa,n,k,s,zarib,pol2,cmpt,t1,t2,b1,b2,JJ,nopv,g,h,a,subb,ff,r,fac:
    a:=[0]:
    t1,b1:=kernelopts(cputime,bytesused):
    nopv:=nops(vars):
     a[1]:=RootOf(J[1]):
     for j from 2 to nopv do
        a:=[op(a),RootOf(subs(seq(x[k]=a[k],k=1..j-1),x[j]=ttt,J[j]),ttt)]:
     od:
     subb:=seq(x[i]=a[i] , i=1..nopv):
    if not issqfr(f) then 
        print("The polynomial is not Squarefree.");
        fff:=subs(subb,f):
        ggg:=evala(Gcd(fff,diff(fff,y))):
        hhh:=evala(Quo(fff,ggg,y)):
        ff:=subs(seq(a[nopv-i+1]=x[nopv-i+1] , i=1..nopv) , hhh);
      else
        ff:=f:  
    fi:
    tord:=plex(y,op(vars)):
    G:=[op(J),ff]:
    flag:=false:cmpt:=0:
    while flag=false do print(cmpt);
        r:=randpoly([op(vars), y], homogeneous, degree = 1);cmpt:=cmpt+1:
        ns, rv := NormalSet(G, tord):
        M := MultiplicationMatrix(r, ns, rv, G, tord);
        chp:=charpoly(M,lambda): 
        if issqfr(subs(lambda=y,chp)) then
            flag:=true:
        fi:
     od:
     fac:=convert(factor(chp),list):
     JJ:=J:
     
     ff:=subs(subb,ff):
     N:=NULL:
     for g in fac do 
        if degree(g,lambda) <>0 then
            h:=subs(lambda=r,subb,g):
            N:=N,evala(Gcd(ff,h)):
        else
            N:=N,g:
        fi:
     od:
     t2,b2:=kernelopts(cputime,bytesused):
     print(t2-t1,b2-b1);
     if {N}={1} or f in {N} then
        RETURN("The polynomial is irreducible");
     fi:
     RETURN(subs(seq(a[nopv-i+1]=x[nopv-i+1] , i=1..nopv) , [N]));
     end:

#################Examples####################
print("******************f1******************");
f1:=expand((y+x[1]*x[3]+x[2]+x[1])*(y-2*x[2]^2+x[3]^2+1)*(y^2+x[1]*x[2]+x[3])):
J:=Basis([x[1]^2-x[2]*x[1]+x[3]*x[1]-x[1]-x[2],x[1]^2-x[3]*x[1]+x[1]-x[2]^2-x[3]*x[2]-x[2]+x[3]^2,-1+x[1]^2+x[3]*x[1]+x[1]-x[2]^2-x[2]+x[3]^2-x[3]],plex(x[3],x[2],x[1]));
A := stan(f1,J, [x[3],x[2], x[1]]);
print("******************f2******************");
f2:=expand((y+x[1])*(y-2*x[4])*(y+x[2]+x[3])):
J:=Basis([x[1]^2+x[3]*x[1]-x[1]*x[4]+x[2]^2-x[2]+x[3]*x[4]-x[4]^2-x[4],x[1]^2+x[2]*x[1]-x[1]*x[4]+x[2]^2+x[3]*x[2]+x[2]+x[4]^2-x[4]
          , 1+x[1]*x[2]-x[1]*x[3]+x[1]*x[4]+x[1]+x[2]^2-x[3]*x[2]+x[2]*x[4]-x[2]+x[3]^2+x[3]*x[4]-x[4]^2 , x[1]^2+x[1]*x[2]+x[3]*x[1]+x[1]*x[4]-x[2]^2+x[3]*x[2]-x[2]+x[3]^2-x[3]*x[4]-x[4]^2+x[4]],plex(x[4],x[3],x[2],x[1])):
A := stan(f2,J, [x[4],x[3],x[2], x[1]]);
print("******************f3******************");
f3:=expand((y-2*x[4]^2+x[3]*x[1]+x[2]+1)*(y+x[2]^2+x[3]*x[4]+x[1]*x[3]+2)):
J:=Basis([-1-x[1]^2+x[3]*x[1]+x[2]^2-x[3]*x[2]+x[3]^2-x[3]*x[4]+x[4]^2+x[4] , 1+x[2]*x[1]+x[3]*x[1]+x[1]-x[3]*x[2]-x[2]*x[4]+x[2]-x[3]^2-x[3]*x[4]-x[3] , 1+x[1]*x[3]+x[1]*x[4]+x[1]+x[2]^2+x[2]*x[4]-x[2]-x[3]-x[4]^2 , 
     x[1]^2+x[1]*x[2]+x[3]*x[1]+x[1]*x[4]-x[1]-x[2]*x[4]-x[3]^2+x[3]*x[4]-x[3]+x[4]^2+x[4]],plex(x[4],x[3],x[2],x[1])):
A := stan(f3,J, [x[4],x[3],x[2], x[1]]);
print("******************f4******************");
f4 := expand((y^2+(x[1]-x[2])*y+x[2])*(y+x[2])*(y-x[4]*x[3])):
J := Basis([x[1]^2+1, x[2]^2+x[1], x[3]^2+x[2], x[4]^2+x[3]], plex( x[4], x[3], x[2], x[1])):
A := stan(f4,J, [x[4],x[3],x[2], x[1]]);
print("******************f5******************");
f5 := expand((y^2+(x[1]-x[2])*y+x[2])^4*(y+x[2])^3*(y-x[4]*x[3])):
J := Basis([x[1]^2+1, x[2]^2+x[1], x[3]^2+x[2]+x[1], x[4]^2+x[3]+x[2]+x[1]], plex( x[4], x[3], x[2], x[1])):
A := stan(f5,J, [x[4],x[3],x[2], x[1]]);
print("******************f6******************");
f6 := expand((y+x[1]*x[3]+x[2]+1)*(y+x[3]*x[4]+x[1]*x[3])*(y^3+x[3]*x[4]+x[1]*x[3])):
J := Basis([x[1]^2+1, x[2]^2+x[1], x[3]^2+x[2], x[4]^2+x[3]+x[2]], plex(x[4], x[3], x[2], x[1])): 
A := stan(f6, J, [x[4], x[3], x[2], x[1]]);