// ################ Introduction // This file contains the Magma 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. // Those parts of this code which are related to compute a // comprehensive Grobner system, is written by D. Kapur, D. Wang and Y. Sun during // their own paper "An efficient algorithm for computing a comprehensive Grobner system of a parametric polynomial system", 2013. // ################ How to execute this file: // 1) Save this file on your system with .txt suffix. // 2) In a Maple worksheet, load this file by: > load "the path you have saved this file/name of // file.txt"; for instance, > load "E:/IntervalGrobSys.txt"; // // 3) Run the main procedure as the following example: // Define the polynomial ring: // P:=PolynomialRing(RationalField(),6,"elim",[1,2],[3,4,5,6]); // where a,b,c,d are parameters which are assigned to intervals. // Define the interval polynomial system: // IPS:=[a*x+b*y+1,b*y^2*x^2+c*x+y+d]; // Define the number of variables and initial null and non-null // conditions as follows: // v:=2; // NC:=[]; // NNC:=[]; // Define also the list of intervals: // Ints:=[[1,2],[-3,3],[2,3],[1,2]]; // Now, run the main procedure: // IntGrobSys(P, v, NC, NNC, IPS, Ints); // For Further help, please contact: Benyamin.M.Alizadeh@gmail.com. Rand:={-10..10}; times:=1; //number of Kapur's check stimes:=times*2; // function NormalSet(P, varlen, GB) //leading monomials lm:=[]; for g in GB do m:=LeadingMonomial(g); Append(~lm, m); end for; //normalset for gb prep:=0; for i:=varlen+1 to Rank(P) do prep:=prep+P.i; end for; prep:=prep+1; prep:=NormalForm(prep, lm); p:=prep; ns:=Monomials(p); newnum:=#ns; repeat oldnum:=newnum; p:=NormalForm(p*prep, lm); ns:=Monomials(p); newnum:=#ns; until oldnum eq newnum; return ns; end function; function Poly_to_Row(p, ns) cp:=Coefficients(p); mp:=Monomials(p); np:=1; res:=[RationalField()|]; for i:=1 to #mp do while ns[np] ne mp[i] do np:=np+1; Append(~res, 0); end while; np:=np+1; Append(~res, cp[i]); end for; while np le #ns do np:=np+1; Append(~res, 0); end while; return res; end function; function MultiplicationMatrix(mr, ns, gb) mtr:=[]; for i:=1 to #ns do nf:=NormalForm(ns[i]*mr, gb); Append(~mtr, Poly_to_Row(nf, ns)); end for; return Matrix(mtr); end function; function Subs(f, v) cf:=Coefficients(f); mf:=Monomials(f); res:=0; for i:=1 to #cf do res:=res+cf[i]*v^(Degree(mf[i])); end for; return res; end function; //only the para part function Zero_Radical(P, varlen, GB) newps:=[]; uni:=[]; for p in GB do flag,_,i:=IsUnivariate(p); if flag then Append(~uni, i); end if; np:=SquarefreePart(p); if np ne p then Append(~newps, np); end if; end for; if #newps ne 0 then newGB:=GroebnerBasis(GB cat newps); else newGB:=GB; end if; for i:=varlen+1 to Rank(P) do if i in uni then continue; end if; ns:=NormalSet(P, varlen, newGB); M:=MultiplicationMatrix(P.i, ns, newGB); chp:=CharacteristicPolynomial(M); nchp:=SquarefreePart(chp); if nchp ne chp then newGB:=GroebnerBasis(newGB cat [Subs(nchp, P.i)]); end if; end for; return newGB; end function; //outputs: 1 for not belonging to gb; 0 for belonging to gb; -1 for zero factors function Radical_test(r, ns, gb) M:=MultiplicationMatrix(r, ns, gb); chp:=CharacteristicPolynomial(M); //not belong to if Coefficient(chp, 0) ne 0 then return 1; end if; t:=LeadingTerm(chp); if t eq chp then //belong to return 0; else //zero factors return -1; end if; end function; //seperate the var-polynomials and para-polynomials procedure SelectPara(P, varlen, PS, ~varPS, ~paraPS) varPS:=[]; paraPS:=[]; for p in PS do if P.varlen gt LeadingMonomial(p) then Append(~paraPS, p); else Append(~varPS, p); end if; end for; end procedure; //check consistent: 4 steps procedure ConsistentCheck(P, varlen, ~ZE, NE, Ints, ~flag, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck) I:=ideal; dim,rrr:=Dimension(I); for i:=varlen+1 to dim do if P.i in ZE then if Ints[i-varlen][1] gt 0 or Ints[i-varlen][2] lt 0 then flag:=false; return; end if; end if; end for; total:=total+1; //trivial-step 1 if #NE eq 0 or #ZE eq 0 then flag:=true; //"trivial check"; tricheck:=tricheck+1; return; end if; s:=1; //the product of all elements in NE for p in NE do s:=NormalForm(s*p, ZE); if s eq 0 then flag:=false; //"trivial check FALSE!!!!!"; tricheck:=tricheck+1; return; end if; end for; //get the dimension HM:=[]; for p in ZE do Append(~HM, LeadingMonomial(p)); end for; I:=ideal

; dim,lst:=Dimension(I); dim:=dim-varlen; lst:=lst[varlen+1..#lst]; //0-dim if dim eq 0 then ZE:=Zero_Radical(P, varlen, ZE); if NormalForm(s, ZE) eq 0 then flag:=false; //"zero-dim check, FALSE!!!!!"; zdcheck:=zdcheck+1; return; else flag:=true; //"zero-dim check"; zdcheck:=zdcheck+1; return; end if; end if; //high-dim-step 2-kapur check for i:=1 to times do //specialization repeat sp:=[]; for j in lst do Append(~sp, P.j - Random(Rand)); end for; new:=NormalForm(s, sp); until new ne 0; spZE:=[]; for p in ZE do newp:=NormalForm(p, sp); if newp ne 0 then Append(~spZE, newp); end if; end for; if #spZE eq 0 then flag:=true; "impossible chack"; continue; end if; //kapur check spZE:=GroebnerBasis(spZE cat sp); if spZE eq [1] then //"bad specialization 1"; continue; end if; num:=0; for p in spZE do if IsUnivariate(LeadingMonomial(p)) then num:=num+1; end if; end for; if num ne Rank(P)-varlen then //"bad specialization 2"; continue; end if; //check whether new belongs to the radical ideal of spZE ns:=NormalSet(P, varlen, spZE); val:=Radical_test(new, ns, spZE); if val ne 0 then flag:=true; kcheck:=kcheck+1; return; end if; end for; //sun check-step 3 new:=s; for i:=1 to stimes //which is the number n of 2^n do res:=0; T:=Terms(new); if #T gt 20 then break; end if; for p in T do res:=res+NormalForm(new*p, ZE); end for; if res eq 0 then flag:=false; //"sun check FALSE!!!!!",i; scheck:=scheck+1; return; end if; new:=res; end for; //damage control - step 4 GB:=GroebnerBasis(ZE cat [P.varlen*s-1]); if GB eq [1] then flag:=false; //"normal check FALSE!!!!!"; ncheck:=ncheck+1; return; else flag:=true; //"normal check TRUE~~~~~"; ncheck:=ncheck+1; end if; end procedure; //make the paraPS radical procedure RadicalOperation_for_paraPS(P, varlen, ~paraPS) newparaPS:=[]; change:=false; for p in paraPS do h:=SquarefreePart(p); if h ne p then change:=true; end if; Append(~newparaPS, h); end for; if not change then return; end if; paraPS:=GroebnerBasis(newparaPS); end procedure; //dim is the upper bound of the dimension & no consistent check procedure RadicalOperation(P, varlen, ~varPS, ~paraPS, ~dim) newparaPS:=[]; change:=false; num:=0; for p in paraPS do h:=SquarefreePart(p); if h ne p then change:=true; end if; Append(~newparaPS, h); if IsUnivariate(LeadingMonomial(h)) then num:=num+1; end if; end for; dim:=Rank(P)-varlen-num; //0-dim if dim eq 0 then if change then newparaPS:=GroebnerBasis(newparaPS); end if; newparaPS:=Zero_Radical(P, varlen, newparaPS); change:=false; for p in newparaPS do if p notin paraPS then change:=true; break; end if; end for; if not change then return; end if; paraPS:=GroebnerBasis(newparaPS); GB:=GroebnerBasis(varPS cat paraPS); SelectPara(P, varlen, GB, ~varPS, ~paraPS); return; end if; //non 0-dim if not change then return; end if; paraPS:=GroebnerBasis(newparaPS); GB:=GroebnerBasis(varPS cat paraPS); SelectPara(P, varlen, GB, ~varPS, ~paraPS); end procedure; //ture for having difference; false for identical. procedure CheckDifference(P, varlen, ZE, newZE, ~flg, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck) total:=total+1; //trival case if #ZE eq 0 then if #newZE eq 0 then tricheck:=tricheck+1; flg:=false; return; else tricheck:=tricheck+1; flg:=true; return; end if; end if; //trival check flag:=true; //ideal identical for p in newZE do if p notin ZE then flag:=false; //ideal not identical break; end if; end for; if flag then //ideal identical tricheck:=tricheck+1; flg:=false; return; //v identical and do not have difference end if; //get the dimension HM:=[]; for p in ZE do Append(~HM, LeadingMonomial(p)); end for; I:=ideal

; dim,lst:=Dimension(I); dim:=dim-varlen; lst:=lst[varlen+1..#lst]; //0-dim if dim eq 0 then zdcheck:=zdcheck+1; flg:=true; return; //ZE is radical ideal end if; HM:=[]; for p in newZE do Append(~HM, LeadingMonomial(p)); end for; I:=ideal

; dim2,_:=Dimension(I); dim2:=dim2-varlen; if dim ne dim2 then tricheck:=tricheck+1; flg:=true; return; //have difference for different dimension end if; //kapur check for polynomials for i:=1 to times do //specialization repeat sp:=[]; for j in lst do Append(~sp, P.j - Random(Rand)); end for; PolyLst:=[]; for p in newZE do new:=NormalForm(p, sp); if new ne 0 then Append(~PolyLst, new); end if; end for; until #PolyLst ne 0; spZE:=[]; for p in ZE do newp:=NormalForm(p, sp); if newp ne 0 then Append(~spZE, newp); end if; end for; if #spZE eq 0 then //"impossible check"; flg:=true; return; end if; //check bad specialization spZE:=GroebnerBasis(spZE cat sp); if spZE eq [1] then //"bad specialization 1"; continue; end if; num:=0; for p in spZE do if IsUnivariate(LeadingMonomial(p)) then num:=num+1; end if; end for; if num ne Rank(P)-varlen then //"bad specialization 2"; continue; end if; //kapur check //check whether all the polynomials in PolyLst belong to the radical ideal of spZE ns:=NormalSet(P, varlen, spZE); for p in PolyLst do val:=Radical_test(p, ns, spZE); if val ne 0 then kcheck:=kcheck+1; flg:=true; return; end if; end for; end for; //sun check-step 3 flag:=true; //assume identical for p in newZE do if stimes eq 0 then break; end if; s:=NormalForm(p, ZE); if s eq 0 then continue; end if; for i:=1 to stimes //which is the number n of 2^n do sum:=0; T:=Terms(s); if #T gt 20 then flag:=false; break; end if; for t in T do sum:=sum+NormalForm(s*t, ZE); end for; if sum eq 0 // p belongs to the radical ideal of ZE then break; end if; s:=sum; end for; if s ne 0 then flag:=false; break; end if; if not flag then break; end if; end for; if flag and stimes ne 0 //all p's in newZE belong to the radical ideal of ZE then scheck:=scheck+1; flg:=false; return; //no difference end if; //damage control - step 4 for p in newZE do s:=NormalForm(p, ZE); if s eq 0 then continue; end if; GB:=GroebnerBasis(ZE cat [P.varlen*s-1]); if GB ne [1] //p do not belong to the radical ideal of ZE then //"normal check True!!!!!!!!"; ncheck:=ncheck+1; flg:=true; return; end if; end for; //"normal check false!!!!!!!!"; ncheck:=ncheck+1; flg:=false; return; //all p's in newZE belong to the radical ideal of ZE //"over"; end procedure; //assume ZE/paraPS is consistent procedure ConsistentCheck_V(P, varlen, ZE, paraPS, NE, Ints, ~flg, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck) I:=ideal; dim,rrr:=Dimension(I); for i:=varlen+1 to dim do if P.i in ZE then if Ints[i-varlen][1] gt 0 or Ints[i-varlen][2] lt 0 then flag:=false; return; end if; end if; end for; total:=total+1; //tirval case if #ZE eq 0 then tricheck:=tricheck+1; flg:=true; return; end if; if #NE eq 0 then tricheck:=tricheck+1; flg:=true; return; //by the assumption end if; //convertion s:=1; //the product of all elements in NE for p in NE do s:=NormalForm(s*p, ZE); end for; if s in BaseRing(P) then tricheck:=tricheck+1; flg:=true; return; //by the assumption end if; newZE:=[]; for p in paraPS do np:=NormalForm(p*s, ZE); if np ne 0 then Append(~newZE, np); end if; end for; if #newZE eq 0 //all polynomials in paraPS be reduced to 0 by ZE then tricheck:=tricheck+1; flg:=false; return; end if; //trival check flag:=true; //ideal identical for p in newZE do if p notin ZE then flag:=false; //ideal not identical break; end if; end for; if flag then //ideal identical tricheck:=tricheck+1; flg:=false; return; //v identical and do not have difference end if; //get the dimension HM:=[]; for p in ZE do Append(~HM, LeadingMonomial(p)); end for; I:=ideal

; dim,lst:=Dimension(I); dim:=dim-varlen; lst:=lst[varlen+1..#lst]; //0-dim if dim eq 0 then zdcheck:=zdcheck+1; flg:=true; return; //ZE is radical ideal end if; //kapur check for polynomials for i:=1 to times do //specialization repeat sp:=[]; for j in lst do Append(~sp, P.j - Random(Rand)); end for; PolyLst:=[]; for p in newZE do new:=NormalForm(p, sp); if new ne 0 then Append(~PolyLst, new); end if; end for; until #PolyLst ne 0; spZE:=[]; for p in ZE do newp:=NormalForm(p, sp); if newp ne 0 then Append(~spZE, newp); end if; end for; if #spZE eq 0 then //"impossible check"; flg:=true; return; end if; //check bad specialization spZE:=GroebnerBasis(spZE cat sp); if spZE eq [1] then //"bad specialization 1"; continue; end if; num:=0; for p in spZE do if IsUnivariate(LeadingMonomial(p)) then num:=num+1; end if; end for; if num ne Rank(P)-varlen then //"bad specialization 2"; continue; end if; //kapur check //check whether all the polynomials in PolyLst belong to the radical ideal of spZE ns:=NormalSet(P, varlen, spZE); for p in PolyLst do val:=Radical_test(p, ns, spZE); if val ne 0 then //"kapur_V check", i; kcheck:=kcheck+1; flg:=true; return; end if; end for; end for; //sun check-step 3 flag:=true; //assume identical for p in newZE do if stimes eq 0 then break; end if; s:=NormalForm(p, ZE); if s eq 0 then continue; end if; for i:=1 to stimes //which is the number n of 2^n do sum:=0; T:=Terms(s); if #T gt 20 then flag:=false; break; end if; for t in T do sum:=sum+NormalForm(s*t, ZE); end for; if sum eq 0 // p belongs to the radical ideal of ZE then break; end if; s:=sum; end for; if s ne 0 then flag:=false; break; end if; if not flag then break; end if; end for; if flag and stimes ne 0 //all p's in newZE belong to the radical ideal of ZE then scheck:=scheck+1; flg:=false; return; //no difference end if; //damage control - step 4 for p in newZE do s:=NormalForm(p, ZE); if s eq 0 then continue; end if; GB:=GroebnerBasis(ZE cat [P.varlen*s-1]); if GB ne [1] //p do not belong to the radical ideal of ZE then //"normal check True!!!!!!!!"; ncheck:=ncheck+1; flg:=true; return; end if; end for; ncheck:=ncheck+1; flg:=false; return; //all p's in newZE belong to the radical ideal of ZE //"over"; end procedure; function HeadCeoff(P, varlen, p, v) g:=p; for i:=1 to varlen do if v[i] ne 0 then g:=Coefficient(g, i, v[i]); end if; end for; return g; end function; //c1 < c2 function Smaller1(c1, c2) if c1[#c1] gt c2[#c2] then return false; else if c1[#c1] lt c2[#c2] then return true; else if #c1 ge #c2 then return true; else return false; end if; end if; end if; end function; //c1 < c2 function Smaller(c1, c2) if #Monomials(c1[#c1]) lt #Monomials(c2[#c2]) then return true; else if #Monomials(c1[#c1]) gt #Monomials(c2[#c2]) then return false; else return c1[#c1] lt c2[#c2]; end if; end if; end function; procedure GetKeyPS(P, varlen, varPS, ~KeyPS, ~NullConditions) KeyPS:=[]; NullConditions:=[]; M:=[]; V:=[]; C:=[]; for i:=1 to #varPS do m:=LeadingMonomial(varPS[i]); v:=Exponents(m); Append(~V, v); CF:=[]; fct:=Factorization(HeadCeoff(P, varlen, varPS[i], v)); for j:=1 to #fct do Append(~CF, fct[j][1]); end for; Append(~C, Sort(CF)); e:=[v[j]:j in [1..varlen]]; for j:=varlen+1 to Rank(P) do Append(~e, 0); end for; Append(~M, Monomial(P, e)); end for; KM:={}; Q:={}; for i:=1 to #M do flag:=false; for m in KM do if IsDivisibleBy(M[i], m) then flag:=true; break; end if; end for; if flag then continue; end if; for j:=i+1 to #M do if IsDivisibleBy(M[i], M[j]) then if M[i] ne M[j] or (M[i] eq M[j] and Smaller(C[j], C[i])) then flag:=true; break; end if; end if; end for; if flag then continue; end if; fct:=C[i]; for j:=1 to #fct do Include(~Q, fct[j]); end for; Include(~KM, M[i]); Append(~KeyPS, varPS[i]); end for; NullConditions:=Sort(SetToSequence(Q)); end procedure; function Reduce_for_NE(NE) newNE:=[]; for p in NE do fct:=Factorization(p); for fc in fct do Include(~newNE, fc[1]); end for; end for; return newNE; end function; //convert NE to ZE in 0-dim ideal procedure NEtoZE(r, ns, gb, ~new, ~total, ~zdcheck) total:=total+1; zdcheck:=zdcheck+1; M:=MultiplicationMatrix(r, ns, gb); chp:=CharacteristicPolynomial(M); chp:=SquarefreePart(chp); //not belong to if Coefficient(chp, 0) ne 0 then new:=1; return; end if; new:=LeadingCoefficient(chp); df:=1; for d:=Degree(chp)-1 to 0 by -1 do c:=Coefficient(chp, d); if c eq 0 then df:=df+1; else for i:=1 to df do new:=NormalForm(new*r, gb); end for; new:=new+c; df:=1; end if; end for; end procedure; //deal with the zero dimension case procedure ZeroDimensionOperation(P, varlen, varPS, paraPS, NE, ~zdres, ~total, ~zdcheck) B:=BaseRing(P); //eliminate NE if #NE eq 0 then gb:=paraPS; else s:=1; //the product of all elements in NE for p in NE do s:=NormalForm(s*p, paraPS); if s eq 0 then zdres:={}; return; end if; end for; if s in B then gb:=paraPS; else ns:=NormalSet(P, varlen, paraPS); NEtoZE(s, ns, paraPS, ~new, ~total, ~zdcheck); if new eq 1 then gb:=paraPS; else gb:=GroebnerBasis(paraPS cat [new]); end if; end if; end if; //pre-operation Done:=[]; Todo:=[]; for p in varPS do np:=NormalForm(p, gb); if np eq 0 then continue; end if; m:=LeadingMonomial(np); v:=Exponents(m); lc:=HeadCeoff(P, varlen, np, v); if lc in B then Append(~Done, np); else Append(~Todo, np); end if; end for; //main loop TD:=[[gb, Done, Todo]]; res:={}; while #TD ne 0 do Br:=TD[#TD]; TD:=TD[1..#TD-1]; gb:=Br[1]; Done:=Br[2]; Todo:=Br[3]; while #Todo ne 0 do p:=Todo[#Todo]; Todo:=Todo[1..#Todo-1]; np:=NormalForm(p, gb); if np eq 0 then continue; end if; m:=LeadingMonomial(np); v:=Exponents(m); lc:=HeadCeoff(P, varlen, np, v); if lc in B then Append(~Done, np); else ns:=NormalSet(P, varlen, gb); NEtoZE(lc, ns, gb, ~new, ~total, ~zdcheck); if new eq 1 then Append(~Done, np); else //add new branch lc=0 Append(~TD, [GroebnerBasis(gb cat [lc]), Done, Append(Todo, np)]); //deal with lc!=0 Append(~Done, np); gb:=GroebnerBasis(gb cat [new]); end if; end if; end while; Include(~res, [gb, [], Done]); end while; zdres:=res; end procedure; procedure PreFullGB(P, varlen, ~ZE, NE) if #ZE eq 0 or #NE eq 0 or #ZE gt 12 then return; end if; oldZE:=ZE; for p in NE do if TotalDegree(p) gt 3 or #Monomials(p) gt 6 then continue; end if; GB:=GroebnerBasis(ZE cat [P.varlen*p-1]); SelectPara(P, varlen, GB, ~varPS, ~ZE); end for; end procedure; procedure PrintBranch(P, F, br) if br[1] ne [] or br[2] ne [] then fprintf F, "\nAll the following conditions should be satisfied:\n"; for p in br[1] do fprintf F, "%o = 0, ", p; end for; if #br[2] ne 0 then CF:=[]; for p in br[2] do Append(~CF, NormalForm(p, br[1])); end for; s:=1; for p in CF do s:=s*p; end for; if s notin BaseRing(P) then CF:=[]; fct:=Factorization(s); for j:=1 to #fct do Append(~CF, fct[j][1]); end for; for p in CF do fprintf F, "%o <> 0, ", p; end for; end if; end if; end if; if #br eq 4 then fprintf F, "\nAt least one of following conditions should be satisfied:\n"; for p in br[3] do fprintf F, "%o <> 0, ", p; end for; end if; fprintf F, "\nThe corresponding GB is:\n%o\n", br[#br]; end procedure; procedure PrintPGB(P, PS, pgb) F:=Open("output.txt","w"); fprintf F, "The polynomial ring is:\n%o\n\n", P; fprintf F, "The system is:\n%o\n\n", PS; fprintf F, "Total number of branches is %o.\n\n", #pgb; i:=0; for br in pgb do i:=i+1; fprintf F, "\nBranch %o:", i; PrintBranch(P, F, br); end for; delete F; end procedure; //the input ZE and NE are assumed to be consistent; procedure PGB(P, varlen, ZE, NE, PS, Ints, ~pgb, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck, trace, F, ~lvl) if trace then fprintf F, "********************PGB function is called (recursive level %o)********************\n\n", lvl; fprintf F, "the inputs of the algorithm PGBMain:\n"; fprintf F, "E is:\n%o\n", ZE; fprintf F, "N is:\n%o\n", NE; fprintf F, "F is:\n%o\n\n", PS; end if; PreFullGB(P, varlen, ~ZE, NE); if trace then fprintf F, "----------Simplifying of E and N----------\nE is:\n%o\n", ZE; fprintf F, "N is:\n%o\n\n", NE; end if; //part-0, compute the full elinimating Grobner basis GB:=GroebnerBasis(PS cat ZE); if trace then fprintf F, "----------Block GB of E cat F is computed----------\nG is:\n\%o\n\n", GB; end if; if GB eq [1] then pgb:={[ZE, NE, [1]]}; if trace then fprintf F, "A trivial branch is obtained"; PrintBranch(P, F, [ZE, NE, [1]]); fprintf F, "\n********************End of PGB function (recursive level %o)********************\n\n", lvl; end if; return; end if; SelectPara(P, varlen, GB, ~varPS, ~paraPS); RadicalOperation(P, varlen, ~varPS, ~paraPS, ~dim); if trace then fprintf F, "----------simplication----------\n"; fprintf F, "polynomials involve variables (G\\G_r):\n%o\n", varPS; fprintf F, "polynomials only involve parameters (G_r):\n%o\n\n", paraPS; end if; //part-1, checking consistent of the new system res:={}; CheckDifference(P, varlen, ZE, paraPS, ~flag, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck); if trace then fprintf F, "----------detect lost parametric space----------\n"; fprintf F, "V(E)-V(G_r) is NOT empty? %o\n\n", flag; end if; if flag then //has difference ConsistentCheck_V(P, varlen, ZE, paraPS, NE,Ints, ~flag, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck); I:=ideal; dim,rrr:=Dimension(I); for i:=varlen+1 to dim do if P.i in ZE then if Ints[i-varlen][1] gt 0 or Ints[i-varlen][2] lt 0 then flag:=false; return; end if; end if; end for; if trace then fprintf F, "V(E)-V(G_r)-V(N) is NOT empty? %o\n\n", flag; end if; if flag then //consistent res:=res join {[ZE, NE, paraPS, [1]]}; ConsistentCheck(P, varlen, ~paraPS, NE, Ints, ~flag, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck); if trace then fprintf F, "A trivial branch is obtained"; PrintBranch(P, F, [ZE, NE, paraPS, [1]]); fprintf F, "\nV(G_r)-V(N) is NOT empty? %o\n\n", flag; end if; if not flag then //not consistent pgb:=res; if trace then fprintf F, "\n********************End of PGB function (recursive level %o)********************\n\n", lvl; end if; return; end if; end if; end if; //part-2, zero dimenstion case if dim eq 0 then res:=res join {[paraPS, NE, varPS]}; if trace then fprintf F, "----------the constraint parametric space is 0-dim----------\n"; fprintf F, "branches with 0-dim parametric constraits:"; fprintf F, "\nbranch:"; PrintBranch(P, F, [paraPS, NE, varPS]); fprintf F, "\n********************End of PGB function (recursive level %o)********************\n\n", lvl; end if; pgb:=res; return; end if; //part-3, high dimension case GetKeyPS(P, varlen, varPS, ~KeyPS, ~NullConditions); if trace then fprintf F, "----------get the noncomparable set----------\n"; fprintf F, "Noncomparable(G\\G_r):\n%o\n", KeyPS; fprintf F, "Factors of leading coefficients of Noncomparable(G\\G_r):\n%o\n\n", NullConditions; end if; newNE:=Reduce_for_NE(NE cat NullConditions); ConsistentCheck(P, varlen, ~paraPS, newNE, Ints, ~flag, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck); if trace then fprintf F, "----------checking whether V(E)-V(N)-V(h) is empty----------\n"; fprintf F, "V(E)-V(N)-V(h) is NOT empty? %o\n\n", flag; end if; if flag then res:=res join {[paraPS, newNE, KeyPS]}; if trace then fprintf F, "A key branch is obtained by our main theorem"; PrintBranch(P, F, [paraPS, newNE, KeyPS]); fprintf F, "\n\n"; end if; end if; newNE:=NE; for i:= 1 to #NullConditions do newZE:=GroebnerBasis(paraPS cat [NullConditions[i]]); RadicalOperation_for_paraPS(P, varlen, ~newZE); ConsistentCheck(P, varlen, ~newZE, newNE, Ints, ~flag, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck); if trace then fprintf F, "----------discussing %o = 0----------\n", NullConditions[i]; fprintf F, "h_1,...,h_%o are:\n%o\n", #NullConditions, NullConditions; if i eq 1 then fprintf F, "V(E, h_%o)-V(N) is NOT empty? %o\n\n", i, flag; else fprintf F, "V(E, h_%o)-V(N)-V(h_1*...*h_%o) is NOT empty? %o\n\n", i, i-1, flag; end if; end if; if flag then lvl:=lvl+1; PGB(P, varlen, newZE, newNE, varPS, ~temppgb, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck, trace, F, ~lvl); lvl:=lvl-1; res:=res join temppgb; end if; Append(~newNE, NullConditions[i]); newNE:=Reduce_for_NE(newNE); end for; pgb:=res; if trace then fprintf F, "********************End of PGB function (recursive level %o)********************\n\n", lvl; end if; return; end procedure; /////////////////////////////////////////////////////////////////////////////// //main function, ZE = zero equations, NE = none-zero equations, //the formation of output is {[ZE, NE, NEI, GB]}, where NEI = none-zero ideal. /////////////////////////////////////////////////////////////////////////////// function IntGrobSys(P, varlen, ZE, NE, PS, Ints) total:=0; tricheck:=0; zdcheck:=0; kcheck:=0; scheck:=0; ncheck:=0; trace:=false; lvl:=0; if #ZE eq 0 then newZE:=[]; else newZE:=GroebnerBasis(ZE); RadicalOperation_for_paraPS(P, varlen, ~newZE); end if; newNE:=Reduce_for_NE(NE); ConsistentCheck(P, varlen, ~newZE, newNE, Ints, ~flag, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck); if flag then PGB(P, varlen, newZE, newNE, PS, Ints, ~IGS, ~total, ~tricheck, ~zdcheck, ~kcheck, ~scheck, ~ncheck, trace, 1, ~lvl); else "Inconsistent input!"; IGS:={}; end if; return IGS; end function;