################ Introduction
# This file contains the Maple code of the IGS algorithm to compute an interval Groebner 
# system for an interval polynomial system, which is stated in our paper entitled:
# "Interval Groebner system and its applications", By:
# Benyamin M.-Alizadeh, Sajjad Rahmany and Abdolali Basiri. 
################ How to execute this file:
# 1) Save this file on your system with .mpl suffix.
# 2) In a Maple worksheet, load this file by: [> read "the path you have saved this file/name of file.mpl";
#   for instance, [> read "E:/IntervalGrobSys.mpl"; 
# 3) Run the main procedure as: [> IntGrobSys(F, Order_p, Order_x, P, Intervals);
#   where, F is the list of the equations in the parametric forms, Order_p (resp. Order_x) is a monomial ordering on
#   parameters (resp. variables), P is the list of parameters and Intervals is the list of intervals.
#   For instance [> IntGrobSys([a*x+b*y+1,b*y^2*x^2+x+y], plex(a,b), tdeg(x,y), [a,b], [[-1,3],[2,4]]);
#
# For Further help, please contact: Benyamin.M.Alizadeh@gmail.com.
###########################################################################
# We need these libraries:
with(LinearAlgebra):
with(combinat):
with(Groebner):
with(PolynomialIdeals):
with(RootFinding):
with(Ore_algebra):
#########################

#########################
RMS := proc (E, f, T) 
local UT,ff,indf, indE, EE, U, ns, rv, Mf, Chf, LM, PE, N, p, test, lm, d, A, n, EEE, R; 
#option trace;
if E=[] then
    RETURN(true);
fi;
indf := indets(f); 
indE := indets(E); 
EE := E; 
if not `subset`(indf, indE) then 
    U := `minus`(indf, indE); 
    EE := [op(E), seq(u-1, u = U)]; 
end if; 
if IsZeroDimensional(`<,>`(op(EE))) then 
    UT := `minus`({op(T)}, indets(EE)); 
    ns, rv := NormalSet(EE, subs(seq(a = NULL, a = UT), T));
    Mf := MultiplicationMatrix(f, ns, rv, EE, T); 
    Chf := CharacteristicPolynomial(Mf, lambda); 
    if LeadingMonomial(Chf, plex(lambda)) <> Chf then 
        RETURN(false); 
    else 
        RETURN(true); 
    end if; 
end if; 
LM := [seq(indets(LeadingMonomial(e, T)), e = E)]; 
PE := `minus`(choose(`union`(indf, indE)), {`union`(indf, indE), {}}); 
N := NULL; 
for p in PE do 
    test := false; 
    for lm in LM while not test do 
        if `subset`(lm, p) then 
            test := true; 
        end if; 
    end do; 
    if not test then 
        N := N, p; 
    end if; 
end do; 
d := max(seq(nops(A), A = [N])); 
N := select(proc (A) options operator, arrow; evalb(nops(A) = d) end proc, [N]); 
#for A in N do 
A:=N[1];
    n := nops(A); 
    EEE := EE; 
    while not IsZeroDimensional(`<,>`(op(EEE))) do 
        R := [seq(RandomTools[Generate](integer(range = -100000 .. 100000)), i = 1 .. n)]; 
        EEE := subs(seq(A[i] = R[i], i = 1 .. n), EE); 
    end do; 
    ff := subs(seq(A[i] = R[i], i = 1 .. n), f); 
    EEE:=Groebner:-Basis(EEE,T);
    UT := `minus`({op(T)}, indets(EEE)); 
    ns, rv := NormalSet(EEE, subs(seq(a = NULL, a = UT), T));
    #ns, rv := NormalSet(EEE, subs(seq(a = NULL, a = A), T)); 
    Mf := MultiplicationMatrix(ff, ns, rv, EEE, subs(seq(a = NULL, a = A), T)); 
    Chf := CharacteristicPolynomial(Mf, lambda); 
    if LeadingMonomial(Chf, plex(lambda)) <> Chf then 
        RETURN(false); 
    else
        RETURN(RadicalMembership(f,<op(E)>));
    end if; 
#end do; 
RETURN(true); 
end proc:

##############################
##########################
# This procedure computes a minimal Dickson basis:
MDBasis:=proc(FF,T1,T2)
local g,N,L,i,flag,F;
#option trace;
F:=FF;
L:=map(u->LeadingMonomial(u,T2),F);
N:=NULL;
for i from 1 to nops(F) do
    flag:=false;
    for g in [N,op(F[i+1..-1])] while not flag do
        if divide(LeadingMonomial(F[i],T2),LeadingMonomial(g,T2)) then 
            flag:=true;
        fi;
    od;
    if not flag then
        N:=N,F[i];
    fi;
od;
RETURN([N]);
end:
##########################
Refine:=proc(E,N,T)
RETURN(E,subs(0=NULL,[seq(NormalForm(n,E,T),n=N)]));
end:

##################################################################

CheckConsistency2:=proc(E,N,L,U,T1)
#option trace;
local flag,i,Test,H,c,l,g,ctm,Y,HH,h;
global Clist;
for g in Clist do
    if g[1]=[E,N] then 
        RETURN(g[2]);
    fi;
od;
if E=[] or N=[] then
    RETURN(true);
fi;

if E=[] then
    RETURN(true);
    if N<>[1] then
        H:=NULL;
        for i from 1 to nops(L) do
            if L[i][2]=infinity then
                H:=H,U[i]-L[i][1]-cat(p,U[i])^2;
            fi;
            if L[i][1]=-infinity then
                H:=H,U[i]-L[i][2]+cat(p,U[i])^2;
            fi;
            if L[i][1]<>-infinity and L[i][2]<>infinity then
                H:=H,U[i]-L[i][1]+(U[i]-L[i][2])*cat(p,U[i])^2;
            fi;
        od;
            ctm:=1;
            for l from 1 to nops(N) do
                ctm:=ctm*w[l]*N[l]-1;
            od;
            Y:=HasRealRoots3([H,ctm]);
            Clist:=[op(Clist),[[E,N],Y]];
            RETURN(Y);
    else
            Clist:=[op(Clist),[[E,N],true]];
            RETURN(true);
    fi;
fi;

if N=[] then
    H:=op(E);
    for i from 1 to nops(L) do
            if L[i][2]=infinity then
                H:=H,U[i]-L[i][1]-cat(p,U[i])^2;
            fi;
            if L[i][1]=-infinity then
                H:=H,U[i]-L[i][2]+cat(p,U[i])^2;
            fi;
            if L[i][1]<>-infinity and L[i][2]<>infinity then
                H:=H,U[i]-L[i][1]+(U[i]-L[i][2])*cat(p,U[i])^2;
            fi;
        od;
        H:=[H];
        HH:=NULL;
        for h in H do
            if indets(E) intersect indets(h)<>{} then
                HH:=HH,h;
            fi;
        od;
        H:=expand([HH]);
            Y:=HasRealRoots3(H);
            Clist:=[op(Clist),[[E,N],Y]];
            RETURN(Y);
fi;


if member(1,E) then
    RETURN(false);
fi;
if {op(E)}={op(N)} then
	RETURN(false);
fi;
Test:=0;
for i from 1 to nops(N) while Test=0 do
    if NormalForm(N[i],E,T1) = 0 then
        Test:=1;
    fi;
od;
if Test=0 then
    flag:=true;
    for i from 1 to nops(N) while flag do 
        if RadicalMembership(N[i],<op(E)>) then
            flag:=false;
        fi;
    od;
    if flag then
        H:=op(E); 
        for i from 1 to nops(L) do
            if L[i][2]=infinity then
                H:=H,U[i]-L[i][1]-cat(p,U[i])^2;
            fi;
            if L[i][1]=-infinity then
                H:=H,U[i]-L[i][2]+cat(p,U[i])^2;
            fi;
            if L[i][1]<>-infinity and L[i][2]<>infinity then
                H:=H,U[i]-L[i][1]+(U[i]-L[i][2])*cat(p,U[i])^2;
            fi;
        od;
        H:=[H];
        HH:=NULL;
        for h in H do
            if indets(E) intersect indets(h)<>{} then
                HH:=HH,h;
            fi;
        od;
        H:=[HH];
        ctm:=0;
        #if N<>[1] then
        #    ctm:=1;
        #    for l from 1 to nops(N) do
        #        ctm:=ctm*w[l]*N[l]-1;
        #    od;
        #    Y:=HasRealRoots3([op(H),ctm]);
        #    Clist:=[op(Clist),[[E,N],Y]];
        #    RETURN(Y);
        #else 
            Y:=HasRealRoots3(H);
            Clist:=[op(Clist),[[E,N],Y]];
            RETURN(Y);
        #fi;
    fi;
fi;
RETURN(false);
end:
#####################################
##########################
# This procedure is the main one:
IntervalPGBmain:=proc(E,N,F,T1,T2,U,L)
local EE,A,PGB,G,NN,g,Gr,Gm,H,h,k,i,G1,T;
#option trace;
PGB:=NULL;
if CheckConsistency2(E,N,L,U,T1) = false then
    RETURN();
fi;
G:=Basis([op(F),op(E)],prod(T2,T1));
if 1 in G then
    RETURN([E,N,{1}]);
fi; 
NN:=NULL;
for g in G do 
    if indets(g) subset {op(U)} then
        NN:=NN,g;
    fi;
od;
Gr:=[NN];
if CheckConsistency2(E,[seq(seq(Gr[i]*N[j],i=1..nops(Gr)),j=1..nops(N))],L,U,T1) then
    if E<>[] or [seq(seq(Gr[i]*N[j],i=1..nops(Gr)),j=1..nops(N))]<>[] then
    	EE,NN:=Refine(E,[seq(seq(Gr[i]*N[j],i=1..nops(Gr)),j=1..nops(N))],T1);
        PGB:=PGB,[EE,NN,{1}];
    fi;
fi;
if CheckConsistency2(Gr,N,L,U,T1)=false then
	RETURN(PGB);
else
    Gm:=MDBasis([op({op(G)} minus {op(Gr)})], T1,T2);
fi;
H:={seq(LeadingCoefficient(g,T2),g=Gm)};
h:=lcm(op(H));
if CheckConsistency2(Gr,[seq(n*h,n=N)],L,U,T1)=true then
	Gr,NN:=Refine(Gr,[seq(n*h,n=N)],T1);
    PGB:=PGB,[Gr,NN,Gm];
fi;
k:=nops(H);
for i from 1 to k do
    if i=1 then
        h:=1;
    else 
        h:=product(H[j],j=1..i-1);
    fi;
    PGB:=PGB,IntervalPGBmain([op(Gr),H[i]],[seq(n*h,n=N)],[op({op(G)} minus {op(Gr)})],T1,T2,U,L);
od;
RETURN(PGB);
end:
###########################
IntGrobSys:=proc(F,T1,T2,U,ListInts)
global Clist;
Clist:=[];
RETURN(IntervalPGBmain([],[1],F,T1,T2,U,ListInts));
end:
########################
########################
with(RegularChains):
with(ParametricSystemTools):
with(SemiAlgebraicSetTools):
HasRealRoots2:=proc(FF)
local R,d,flag,B,i,infolevel,F;
#option trace;
d:=nops(indets(FF));
F:=convert(expand(subs(seq(w[i]=cat(w,i),i=1..d),FF)),rational);
infolevel[RegularChains] := 1:
R := PolynomialRing([op(indets(F))]);

flag:=false;
for i from 1 to d-1 while not flag do 
    B:=RealRootClassification(F,[],[],[],i,1..k,R):
    if nops(B[1])<>0 then
        flag:=true;
    fi;
od;
RETURN(flag);
end:
##########################
with(RAGMaple):
with(RootFinding):
HasRealRoots3:=proc(FF)
local f,F,HRS;
#option trace;
#RETURN(HasRealRoots2(FF));
#RETURN(HasRealRoots(FF));
#F:=[seq(f=0,f=Basis(FF,plex(op(indets(FF)))))];
F:=[seq(f=0,f=FF)];
HRS:=HasRealSolutions(F);
if HRS=[] then
    RETURN(false);
fi;
RETURN(true);
end: