# This file contains the Maple code of three algorithms: 
# The first is DetNorPos to transform a zero-dimensional
# radical ideal into normal position.
# The second is "ImpDetNorPos" which is an improvement of the first. 
# Finally, the third is "ImpDetNorPosNR" which is devoted to transform
# a zero-dimensional non-radical ideal into normal position.
# These algorithms are all presented in a paper authored by: 
# Benyamin M.-Alizadeh and Amir Hashemi, 2019. 
# If there is any comment, please contact with Benyamin.M.Alizadeh@gmail.com.
################################
with(PolynomialTools):with(LinearAlgebra):
with(PolynomialIdeals):with(Groebner):with(RandomTools):
################################
# Deterministic Normal Position : The Maple code of Algorithm 3
################################
DetNorPos := proc (F, var) 
#option trace;
local nV, nslist, f, V, Tordd, G, ns, rv, flag, n, Flag, i, W,t1,t2,b1,b2; 
if IsRadical(<op(F)>) = false then 
   error "For this version, the input ideal must be radical. Please use ImpDetNorPosNR.";
end if; 
t1,b1:=kernelopts(cputime,bytesused);
f := 0; 
V := indets(F) minus {var};
nV:=nops(V); 
Tordd := plex(seq(V[nV-i+1],i=1..nV), var); 
G := Basis(F, Tordd); 

#G:=Invg2v(F,Tordd);
ns, rv := NormalSet(G, Tordd); 
nslist:=convert(ns,list);
ns:=sort(nslist,(p,q)->TestOrder(q,p,Tordd)); 
flag := false; 
n := nops(ns); 
while not flag do 
    Flag := false; 
    for i from 1 to n-1 while not Flag do 
        W := indets(ns[i]); 
        if W <> {var} then 
            f := f+(W minus {var})[-1]; 
            #G := Basis(subs(var = var+(W minus {var})[-1], G), Tordd); 
            G:=Invg2v(expand(subs(var = var+(W minus {var})[-1], G)), Tordd);
            ns, rv := NormalSet(G, Tordd); 
            nslist:=convert(ns,list);
            ns:=sort(nslist,(p,q)->TestOrder(q,p,Tordd));
            Flag := true; 
        end if; 
    end do; 
    if not Flag then 
        flag := true; 
    end if; 
end do;
t2,b2:=kernelopts(cputime,bytesused); 
printf("%-1s %1s %1s %1s          : %1a %3a\n",The, cpu, time, is,t2-t1,(sec)):
 printf("%-1s %1s %1s          : %5a %3a\n",The,used,memory,b2-b1,(bytes)):
 printf("%-1s %1s %1s           : %1a\n",Cardinal,of,NS,n):
 printf("%-1s %1s %1s        : %2a\n",The,linear,change,var=var+f):
  printf("%-1s %1s %1s    : %2a\n",Correctness,of,result,evalb(degree(G[1])=n)):
RETURN(G); 
end proc:
#####################
# Improved Version : The Maple code of Algorithm 4
#####################
IsDep:=proc(F)
global G, n, Ns,Tordd,S;
local f, H;
#option trace;
f:=add(a[i]*F[i],i=1..nops(F));
f:=subs(S,f);
H:=[seq(coeff(f,y[i]),i=1..n-1),subs(seq(y[i]=0,i=1..n-1),f)];
H:=Basis(H,tdeg(seq(a[i],i=1..nops(F))));
if {seq(h/LeadingCoefficient(h,tdeg(seq(a[i],i=1..nops(F)))),h=H)} = {seq(a[i],i=1..nops(F))} then
#if solve(H) = {seq(u=0,u=indets(H))} then
    RETURN(false);
else
    RETURN(true);
fi;
end:
#####################
ImpDetNorPos := proc (F, var) 
#option trace;
global G,Ns,Tordd,n,S;
local f, V, rv, flag, ns, Flag, i, W, t1, t2, b1, b2, nV, K, temp, m, cmpt, u;
if IsRadical(<op(F)>) = false then 
    error "For this version, the input ideal must be radical. Please use ImpDetNorPosNR."; 
end if; 
f := 0; 
V := indets(F) minus {var};
nV:=nops(V); 
Tordd := plex(seq(V[nV-i+1],i=1..nV), var); 
G := Basis(F, Tordd); 
t1,b1:=kernelopts(cputime,bytesused);
ns, rv := NormalSet(G, Tordd); 
n := nops(ns); 
Ns:=sort(ns,(p,q)-> TestOrder(q,p,Tordd));
S:=seq(Ns[i]=y[i],i=1..n-1);
Flag:=false;
while not Flag do 
    flag:=false;
    K:=[1];
    temp:=1;
    for i from 1 to n-1 while not flag do 
        temp:=NormalForm(temp*(var-f),G,Tordd);
        K:= [op(K),temp];
        if IsDep(K) then
            flag:=true;
        fi;
    od;
    m:=i-2;
    if not flag or m=n then
        Flag:=true;
    else 
        cmpt:=0; 
        while cmpt=0 and V<>{} do
            u:=V[-1];
            if IsDep([op(K[1..-2]),NormalForm(u,G,Tordd)]) then
                V:=V[1..-2];
            else
                f:=f+u; 
                cmpt:=1;
            fi;
        od;
    fi;
od;
t2,b2:=kernelopts(cputime,bytesused); 
printf("%-1s %1s %1s %1s          : %1a %3a\n",The, cpu, time, is,t2-t1,(sec)):
 printf("%-1s %1s %1s          : %5a %3a\n",The,used,memory,b2-b1,(bytes)):
 printf("%-1s %1s %1s       : %1a\n",Dimension,of,ideal,n):
 printf("%-1s %1s %1s        : %2a\n",The,linear,change,var=var+f):
end proc:
#######################################
# Normal Position for non-radical ideals : The Maple code of Algorithm 5
#######################################
ImpDetNorPosNR := proc (F, var) 
#option trace;
#global G,Ns,Tordd,n,S;
local f, V, rv, flag, ns, Flag, i, W, t1, t2, b1, b2, nV, K, temp, m, cmpt, u, fgold, Tordd, G, n, Ns, dtemp, l, H, p;
if IsRadical(<op(F)>) then 
    print("For radical ideals it is better to use ImpDetNorPos."); 
end if; 

fgold := 0; 
V := indets(F) minus {var};
nV:=nops(V); 
Tordd := plex(seq(V[nV-i+1],i=1..nV), var); 
G := Basis(F, Tordd); 
t1,b1:=kernelopts(cputime,bytesused);
ns, rv := NormalSet(G, Tordd); 
n := nops(ns); 
Ns:=sort(ns,(p,q)-> TestOrder(q,p,Tordd));
for u in V do
    f:=fgold;
    dtemp:=0;
    for l from 1 to 2*n-2 do
        f:=f+u; 
        H:=Basis(subs(var=var+f,G),Tordd);
        p:=SquareFreePart(H[1]);
        if dtemp<degree(p) then
            fgold:=f;
            dtemp:=degree(p);
        fi;
    od;
od;
t2,b2:=kernelopts(cputime,bytesused); 
printf("%-1s %1s %1s %1s          : %1a %3a\n",The, cpu, time, is,t2-t1,(sec)):
 printf("%-1s %1s %1s          : %5a %3a\n",The,used,memory,b2-b1,(bytes)):
 printf("%-1s %1s %1s       : %1a\n",Dimension,of,ideal,n):
 printf("%-1s %1s %1s        : %2a\n",The,linear,change,var=var+fgold):
end proc:

######### Example
F:=[x^4+2, y^2+2*z+2, (z+y)^2+1];
ImpDetNorPos(F, z);
print("#################");
F := [x[6]^3-16*x[6]^2+81*x[6]-126, x[5]-5, (x[4]^2-9*x[4]+20)^2];
ImpDetNorPosNR(F, x[5]);